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Chapter  1 


Introduction 

The  magnitude  and  phase  of  the  Fourier  transform  of  an  arbitrary  multidimensional 
signal  are  independent  functions  of  frequency.  In  many  applications,  however,  there 
is  additional  information  regarding  the  signal  which  provides  a  very  strong  connection 
between  its  Fourier  transform  magnitude  and  phase.  One  example  of  such  additional 
information  is  the  common  condition  that  the  signal  is  non-zero  only  over  a  specified 
region.  In  this  case,  it  has  been  shown  that  almost  all  multidimensional  signals  which 
are  non-zero  only  over  a  specified  region  are  uniquely  specified,  in  a  sense,  by  knowl¬ 
edge  of  only  its  Fourier  transform  magnitude  [1,2).  Hence,  once  the  Fourier  transform 
magnitude  is  known,  the  Fourier  transform  phase  is  determined  as  well.  For  this  rea¬ 
son,  the  problem  of  reconstruction  from  Fourier  transform  magnitude  is  also  called  the 
phase  retrieval  problem. 

The  reconstruction  of  a  two-dimensional  signal  from  its  Fourier  transform  magni¬ 
tude  has  been  the  object  of  much  study.  This  interest  is  guided  by  the  wide  range  of 
applications  of  results  in  this  area.  One  such  application  is  in  astronomy  [3j.  The  effect 


. . . . .  ^  t . . . 


S. 

.  a.  r  . 


7 


limiting  the  resolution  capabilities  of  the  largest  optical  telescopes  is  not  the  diffraction 
limit  of  the  lens,  but  rather  the  turbulence  of  the  earth’s  atmosphere.  During  a  short 
time  period,  the  atmosphere  can  be  considered  to  introduce  on  the  incoming  optical 
wave  a  spatially  varying,  time-invariant  random  phase  delay  due  to  inhomogeneities 
induced  by  thermal  gradients.  These  inhomogeneities  are  slowly  varying  with  respect 
to  short  exposure  times.  Thus,  over  such  a  short  time  period,  the  atmosphere  can  be 
modeled  as  a  glass  plate  of  spatially  varying  thickness  over  the  telescope  aperture. 

This  phase  aberration,  although  it  blurs  each  individual  exposed  image,  does  not 
affect  the  spatial  autocorrelation  function.  A  way  of  circumventing  this  blurring  effect 
is  to  first  measure  an  accurate  estimate  of  the  spatial  autocorrelation  function.  This 
can  be  done  via  Labeyrie  interferometry  [4],  In  this  procedure  an  interferometer  is 
used  to  image  the  spatial  autocorrelation  function  over  a  small  time  period.  This 
estimate  will  be  very  noisy  because  of  the  short  exposure  time.  The  signal-to-noise 
ratio  however  can  be  increased  by  averaging  several  short  exposures.  Thus  a  diffraction 
limited  autocorrelation  function  can  be  measured  which  is  net  affected  by  atmospheric 
blurring.  It  is  clear  that  a  reliable  method  for  extracting  the  image  of  the  astronomical 
object  from  such  interferometer  data  would  in  effect  greatly  increase  the  resolution 
capabilities  of  earth-based  telescopes. 

A  possible  application  of  phase  retrieval  to  electron  microscopy  lies  in  the  possibility 
of  indirect  phase  measurement  from  magnitude  measurement  [5j.  Photographic  film 
can  only  record  the  intensity  of  the  field  impinging  on  it.  However,  the  phase  of  the 
field  provides  important  information  on  the  object  being  viewed.  For  example,  thin 
objects  may  be  considered  as  modulating  the  phase  of  the  electron  wave  while  not 


affecting  the  magnitude.  The  actual  phase  delay  introduced  depends  on  the  thickness 
and  composition  of  the  specimen.  Retrieving  the  phase  delays  from  the  recorded  field 
intensity  would  yield  an  indirect  way  of  measuring  the  specimen  properties. 

X-ray  crystallography  is  a  third  realm  where  reconstruction  from  Fourier  transform 
magnitude  may  prove  useful  [6].  Physical  arguments  show  that  the  angles  at  which  the 
x-rays  are  diffracted  from  a  crystal  specimen  and  the  intensity  of  the  diffracted  wave 
at  each  angle  are  related  to  the  Fourier  transform  magnitude  of  the  electron  density  of 
the  crystal  under  study.  An  important  part  of  crystallography  is  the  task  of  deducing 
the  arrangement  of  atoms  in  the  crystal  from  knowledge  of  such  diffraction  data. 

The  importance  of  the  phase  retrieval  problem  has  led  several  researchers  to  propose 
algorithms  for  reconstruction  from  Fourier  transform  magnitude.  However,  previously 
presented  algorithms  fall  into  either  of  two  categories;  they  are  heuristic  algorithms 
which  often  do  not  converge  to  the  true  reconstruction,  or  they  are  computationally 
too  expensive  for  even  moderate  size  signals.  The  purpose  of  this  thesis  is  to  present 
a  new  algorithm  for  reconstruction  of  multidimensional  discrete  signals  from  Fourier 
transform  magnitude  which  is  a  ciosed  form  solution  to  the  problem  and  which  has 
been  used  with  success  in  reconstructing  signals  of  moderate  size. 

The  thesis  is  divided  into  eight  chapters;  Chapter  2  develops  an  appropriate  nota¬ 
tion,  reviews  basic  properties  of  signals  and  introduces  some  mathematical  concepts. 
Chapter  3  is  concerned  with  the  review  of  previous  results  in  reconstruction  of  both 
one-  and  two-dimensional  signals  from  the  Fourier  transform  magnitude.  The  formu¬ 
lation  of  the  phase  retrieval  probiem  as  a  bivariate  poiynomiai  factorization  problem 
is  given  in  Chapter  4.  Previous  algorithms  for  factoring  polynomials  in  two  variables 
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and  their  connection  and  possible  application  to  the  phase  retrieval  problem  are  also 
discussed.  Based  on  this  framework,  a  new  algorithm  for  factoring  large  polynomials 
in  two  variables  is  developed  in  Chapter  5.  This  leads  to  a  new  closed  form  algorithm 
for  solving  the  phase  retrieval  problem.  Examples  of  the  application  of  the  new  phase 
retrieval  algorithm  is  the  subject  of  Chapter  6.  Chapter  7  discusses  the  application  of 
the  ideas  developed  in  this  thesis  to  the  areas  of  general  bivariate  polynomial  factor¬ 
ization  and  filter  stability  testing.  A  summary  of  the  work  presented  and  suggestions 
for  future  research  is  the  subject  of  Chapter  8. 
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Chapter  2 


Notation  and  Basic  Properties 


In  this  chapter  we  review  some  concepts  from  digital  signal  processing  and  polyno¬ 
mials  in  one  or  several  variables.  The  purpose  of  this  review  is  primarily  to  develop  a 
consistent  notation  and  provide  a  base  for  our  later  development.  The  interested  reader 
may  consult  [7]  for  a  more  detailed  discussion  of  the  signal  processing  topics  described 
here.  A  development  of  the  properties  of  polynomials  reviewed  here  is  contained  in  [8]. 


2.1  One-  and  Two-Dimensional  Signals 

A  one-dimensional  discrete  signal  xjn]  or  a  two-dimensional  discrete  signal  z[m,  n] 
is  a  real  or  complex  function  of  a  single  integer  index  n  or  two  integer  indices  m  and  n 
respectively.  The  support  of  x[n ’  is  the  set  (or  sometimes  a  superset)  of  all  indices  n 
such  that  xfn;  is  non-zero.  The  support  of  x[m,  nj  is  also  defined  as  the  set  of  ail  index 
pairs  such  that  x[m,  n\  is  non-zero. 

In  our  discussion  we  will  find  special  cases  of  support  to  be  especially  useful.  A 


signal,  either  one-  or  two-dimensional,  is  said  to  be  of  finite  extent  if  its  support  is 
bounded;  in  one  dimension,  this  means  that  there  is  an  index  pair,  (nmtn,nroaa)  such 
that  x[n]  =  0  for  n  >  nTO0J  and  n  <  nm,„.  Similarly  in  two  dimensions  a  signal  is  of  finite 
extent  if  there  is  an  index  quadruple  (mm0,,  mmin ,  nmax,  nm,„)  such  that  x[m,  n]  =  0 
whenever  any  of  the  four  conditions  below  hold: 

m  >  ,  77i  <  mmin,  ti  >  n,m0J,  n  <  7imin  (2-1) 

Pictorially,  this  means  that  the  non-zero  values  of  x[m,  nj  can  be  enclosed  in  a  box, 
Figure  2.1. 


n 

* 


Figure  2.1:  A  two-dimensional  signal  with  support  x  [iVmtn, 

We  will  denote  a  region  of  support  consisting  of  all  indices  a  <  n  <  6  as  |a,  6’. 
In  two  dimensions  a  support  comprising  all  index  pairs  (m,  n)  satisfying  a  <  n  <  b, 
c  <  m  <  d  will  be  denoted  by  [a, fcj  x  [c,d\.  We  will  abbreviate  [a,6j  x  [a, 6j  by  [a,  6I-. 
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Associated  with  any  signal  x{n]  is  a  Laurent  series  called  its  z-transform, 

X{z)  =  Y,x\n\z~n  (2.2) 

n 

where  z  is  a  complex  number.  The  set  of  all  z  where  the  summation  converges  is  called 
the  region  of  convergence  (ROC)  of  X(z).  In  this  thesis,  we  will  be  dealing  primarily 
with  signals  of  finite  extent,  in  which  case  the  ROC  includes  all  of  the  complex  z-plane 
with  the  possible  exclusion  of  the  origin  or  infinity.  Evaluating  the  z-transform  in  (2.2) 
on  the  unit  circle  |zj  =  1  yields  the  Fourier  transform, 

AV")  =  £*[n]«"y,M’  (2-3) 

n 

Even  if  x[n]  is  real,  its  Fourier  transform  may  be  complex-valued,  and  can  thus  be 
represented  in  terms  of  its  real  and  imaginary  components, 

X(e>')  =  Xr(ei')  +  jXi(e>')  (2.4) 

or  in  polar  form 

X(e'v)  =  j  A’(e;t')|e;e,(r*  (2.5) 

Two-dimensional  discrete  signals  also  have  a  z-transform  which  is  a  complex-valued 
function  of  two  complex  numbers  w  and  z, 

X(w,  z)  =  Y,  H  x[m>  n\w-mz~n  (2.6) 

m  n 

The  definition  of  an  ROC  also  applies  to  two-dimensional  z-transforms.  Whenever  the 
ROC  includes  the  bicircle  ju;|  =  1,  jzj  =  1,  the  Fourier  transform  of  zfm,  nj  is  defined: 


An  important  signal  which  is  generated  from  z[n|  is  its  autocorrelation  function, 
defined  via  the  summation  below, 

r(nl  =  53  *[/]**{/  +  nj  (2.8) 

; 

It  is  straightforward  to  show  that  if  x[n)  is  a  finite  extent  signal  of  support  [a,  6]  then 
r[n]  has  support  [—(b  —  a),  (b  —  a)).  Note  that  the  autocorrelation  function  is  invariant 
to  multiplication  of  x[n]  by  a  constant  of  unit  magnitude  or  to  a  translation  of  xjnj. 

From  the  convolution  theorem  [7],  the  following  relationship  is  derived  between 
X(z)  and  the  z-transform  of  r[n],  R(z), 

m  =  xw(±)  (2.9) 

When  evaluated  on  the  unit  circle,  the  above  equation  collapses  to  a  relationship  be¬ 
tween  the  Fourier  transforms  of  x[n\  and  r{n], 

R{e>v)  =  \X[e’r)\2  (2.10) 

From  (2.10)  one  can  make  an  important  observation  that  if  two  signals  have  the  same 
autocorrelation  function,  they  must  have  the  same  Fourier  transform  magnitude,  and 
vice  versa. 

The  autocorrelation  function  of  x[m,  nj  is  similarly  defined  by 

r[m,  n]  =  z[fc,/]x*[fc  +  m,  /  +  n]  (2.11) 

i  i 

The  support  of  rim,  nj  is  given  by  !-(6  —  a),  (6  -  a)]  x  [~(d  -  c),(d  -  c)]  for  xjm,  n] 
with  support  [a,b]  x  \c,d}. 
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The  convolution  theorem  also  applies  to  the  two-dimensional  case;  the  z-transform 
of  r[m,n],  R(w,z),  is  given  by 

R(w,z)  =  X(tv,z)X-(±  1)  (2.12) 

z 

On  the  unit  bicircle,  (2.12)  becomes 

R(eJU,eJ')  =  \X[e’u,e’v)\2  (2.13) 

Thus  we  see  that  as  in  the  one-dimensional  case,  if  the  Fourier  transform  magnitude  of 
a  signal  is  known  then  its  autocorrelation  function  is  known  also  and  vice  versa. 


2.2  Polynomials  in  One  and  Two  Variables 

In  the  subsequent  discussion,  we  will  be  dealing  primarily  with  signals  of  finite  ex¬ 
tent.  In  this  case,  the  corresponding  z-transforms  are  essentially  polynomials.  There¬ 
fore,  it  is  important  to  understand  some  properties  of  polynomials  which  will  be  used 
later  on. 

A  polynomial  in  one  variable  z  is  a  function  of  the  form 

.v 

P{z)  =  Y,P"Zn  (2.14) 

n=0 

where  the  p„  are  complex  numbers  and  N  is  finite.  The  degree  of  a  polynomial  in  one 
variable  (or  one-dimensional  polynomials)  is  the  largest  power  to  which  the  indeter¬ 
minate  variable  is  raised.  The  degree  of  a  polynomial  p(z)  will  be  denoted  by  deg(p). 
Thus  the  polynomial  p(z)  in  (2.14)  above  has  degree  deg(p)  =  N 

As  a  result  of  the  Fundamental  Theorem  of  Algebra  [8],  all  one  dimensional  poly¬ 


nomials  of  degree  N  can  always  be  expressed  as  a  product  of  N  polynomials  of  degree 
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1  and  a  constant, 


.v 

P(*)  =  P.Y  IK*”2")  (2.15) 

n=l 

This  product  representation  is  unique  up  to  a  permutation  of  the  product  terms.  Since 
p(zn)  =  0  for  all  n,  z„  is  called  a  zero  of  p{z).  We  note  that  z„  will  generally  be  complex, 
even  if  the  p„  set  is  real. 

A  polynomial  p(z)  with  deg(p)  >  0  is  called  reducible  if  it  can  be  expressed  as  the 
product  of  two  polynomials  pi(z),  p2(z)  with  deg(pi)  >  0  and  deg(p2)  >  0,  i.e., 

p{z)  =  pl{z)p2{z)  (2.16) 

If  no  such  decomposition  is  possible,  then  p(z)  is  called  irreducible.  From  the  decompo¬ 
sition  presented  in  (2.15),  we  see  that  the  only  irreducible  polynomials  in  one  variable 
are  polynomials  of  degree  1  '.  We  will  see,  however,  that  the  situation  is  quite  different 
for  polynomials  in  two  (or  more)  variables. 

Associated  with  any  polynomial  is  a  “mirror”  polynomial  consisting  of  coefficients 
in  reversed  order  and  conjugated.  For  example,  for  the  polynomial  p[z)  in  (2.14),  the 
mirror  polynomial  p(z)  is  defined  by 

p(*)  =  £p.V- (2-17) 

n=0 

There  is  a  very  simple  relationship  between  the  zeros  of  p(z)  and  p(z);  namely,  if  z0  is 
a  zero  of  p(z),  then  z,j-1  is  a  zero  of  p(z). 

•Strictly  speaking,  irredncihility  is  defined  with  respect  to  a  specific  field.  Whether  a  polynomial  is 
irreducible  or  not  may  depend  on  the  field  of  interest.  However,  in  our  discussion  we  will  only  be  dealing 
with  polynomials  over  the  field  of  complex  numbers. 


16 


One  can  also  study  polynomials  in  two  variables,  also  called  two-dimensional  poly¬ 


nomials  or  bivariate  polynomials.  In  this  case  there  are  two  variables  w,  z, 

M  -V 

PK*)  =  II  pm.nwmzn  (2.18) 

m=0  r»=0 

where  again  pm.„  are  allowed  to  be  complex  numbers.  The  degree  in  w,  degu.(p),  of 
p(w,  z)  is  the  highest  order  to  which  the  indeterminate  w  is  raised.  In  the  example 
above,  degir{p)  =  M.  Similarly,  the  degree  in  z  of  p(w,z),  deg:(p),  is  N.  The  degree  of 
the  polynomial  p(w,  z)  deg(p),  is  defined  by  the  pair  of  integers  (degv{p),deg:{p))-,  in 
this  case,  deg{p)  =  ( M,N ).  We  will  also  define  the  total  degree  of  p(u>,  z),  totdeg(p), 
as  the  degree  of  the  univariate  polynomial  p(u>,  tu).  We  note  that  in  some  areas  of 
mathematics,  e.g.,  algebraic  geometry,  the  total  degree  is  considered  the  degree  of 
p{w,z). 

A  decomposition  of  an  arbitrary  bivariate  polynomial  into  a  product  of  polynomials 
of  a  lower  degree  is  not  always  possible.  This  is  because,  unlike  the  case  for  one¬ 
dimensional  polynomials,  there  are  irreducible  polynomials  in  two  variables  for  any 
degree,  except  of  course  polynomials  which  are  of  degree  0  in  one  of  the  variables. 

Example:  The  polynomial  of  degree  (M.N)  below  is  easily  checked  fo  be  irreducible  for 
any  .V  >  0  and  M  >  0  [9], 

p(tmz)  =  p  o  +  piw  +  ■  ■  ■  +  pi«ru  +  (2.19) 

If  a  product  decomposition  does  exist,  then  it  is  essentially  unique  as  expressed  in 
the  following  theorem  10 ', 

Theorem  2.1  Every  polynomial  f(w,z)  is  expressible  as  the  product 


f  =  PiPz---Pk 


(2.20) 


of  a  finite  number  of  irreducible  polynomial  factors  Pi,  and  every  other  such  expression 
for  f  has  factors  which  are  the  same  except  possibly  for  constant  multipliers. 

Note  that  this  theorem  is  equivalent  to  the  Fundamental  Theorem  of  Algebra  except 
that  in  this  case  there  is  no  specification  of  the  degrees  of  each  irreducible  factor  P,. 

An  important  representation  of  a  bivariate  polynomial  p(tt> ,  z)  is  to  consider  it  as  a 
polynomial  in  w  which  has  coefficients  consisting  of  polynomials  in  z, 

p(w,  z)  =  p0(z)  +  pi{z)w  +  •  •  •  pxt(z)ur'f  (2.21) 

If  deg(p)  =  (A/,  N)  then  the  degree  of  each  p„(z)  must  be  less  than  or  equal  to  N.  Of 
course,  one  can  similarly  consider  the  same  polynomial  p(u>,  2)  as  a  one-dimensional 
polynomial  in  2  with  coefficients  consisting  of  polynomials  in  w. 

Similarly,  p(w,z),  the  “mirror”  image  of  p(w,z),  is  defined  by 

~p{w,z)  =  £  £  p.v_m..v_ntnn,2n  (2.22) 

m=  0  n=0 

The  zeros  of  p(u>,  z)  and  p(u>,  z)  are  related  via  a  relationship  analogous  to  the  univariate 
case;  if  (u>0,2o)  is  a  zero  of  p(w,z),  then  (u>q-1,2o-1)  is  a  zero  of  p(w,z). 


2.3  Z-Transforms  of  Finite  Extent  Signals 


Recall  that  a  one-dimensional  discrete  signal  zjn]  with  finite  extent  support  [a,  6] 
has  z-transform, 

b 

X{z)  =  YL  (2-23) 

n=o 

The  expression  above  can  be  written  equivalently  as 

b—n 

X(z)  =  z~tYlz\t>-n}zn  (2.24) 


n=0 
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The  second  term  can  be  identified  as  a  polynomial  in  the  variable  z  with  coefficients 
p„  =  x[b  -  n].  We  will  call  the  expression 

P»(*)  =  E  x\b  ~  nlz"  (2.25) 

n=0 

the  polynomial  associated  with  x[n].  Thus,  the  two  signals  z[n]  and  z[n  +  A;]  for  any 
fixed  k  have  the  same  associated  polynomial.  We  will  assume  that  the  degree  of  p,(z)  is 
as  small  as  possible;  thus,  any  associated  polynomial  of  degree  n  has  a0  #  0  and  a„  ^  0. 
Formally,  we  define  Pi(z)  as  the  polynomial  of  least  degree  such  that  X(z )  =  zkpI(z) 
for  some  integer  k. 

There  is  a  relationship  concerning  the  associated  polynomials  of  x\n ]  and  r[n],  which 
is  very  similar  to  the  convolution  theorem  relationship  of  (2.9), 

Pr(z)  =  P,(z)p*(z)  (2-26) 

where  we  recall  that  p.,(z)  is  the  mirror  polynomial  of  p,(z). 

We  can  also  define  associated  polynomials  for  two-dimensional  signals.  Consider 
the  z-transform  of  a  two-dimensional  signal  x[m,  nj  of  support  [a,  6]  x  [c,  d]t 

X(w,  z )  =  w~dz~h  ^  ]T]  x[d  —  m,b  —  n\wmzn  (2.27) 

m=0  n=0 

Again,  we  define  the  polynomial 

d—t  b—a 

p,(u;,  z)  =  r.  y.  z[d  -  m,  6  -  n\wmzn  (2.28) 

m=0  n=0 

as  the  polynomial  associated  with  xlm ,  nl.  The  corresponding  relationship  between  the 
polynomial  associated  with  xlm,  n]  and  r[m,  n)  is  given  by 

pr(w,  z)  =  p,(w,  z)p,(w,  z )  (2.29) 
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2.4  The  Phase  Retrieval  Problem 


Using  the  notation  presented  in  this  chapter,  we  can  state  succinctly  what  is  the 
specific  problem  we  are  trying  to  solve.  Consider  an  unknown  two-dimensional  discrete 
signal  x[m,  n].  The  problem  is  to  reconstruct  ar[m,  n]  given, 

•  The  support  of  x [m,  n]. 

•  The  Fourier  transform  magnitude,  \X(ein ,eiv)\>  of  x[m,  n],  at  all  frequencies  2. 

This  problem  can  be  reformulated  in  several  different  ways.  From  the  relationship 
of  (2.13),  we  know  that  knowledge  of  the  Fourier  transform  magnitude  is  equivalent  to 
knowledge  of  the  autocorrelation  function  of  r[m,  nj.  Thus  the  phase  retrieval  problem 
can  be  stated  as,  reconstruct  x\ m,  n)  given, 

•  The  support  of  i[m,  n). 

•  The  autocorrelation  function,  r[m,  n],  of  x[m,n}. 

The  polynomial,  p,(u>,  z)  associated  with  i [m,  n]  is  an  almost  invertible  represen¬ 
tation  of  a  signal,  since  from  the  coefficients  of  pz(w,z)  we  can  restore,  to  a  linear 
shift,  the  original  signal.  We  will  find  that  this  shift  ambiguity  is  inherent  in  phase 
retrieval.  Using  associated  polynomials,  the  phase  retrieval  problem  can  be  considered 
as  reconstructing  pz{w,  z )  given, 

•  the  degree  of  pz(w,z). 

JIt  can  he  shown  that  for  finite  extent  signals  it  is  only  required  that  the  Fourier  transform  magnitude 
he  known  on  a  sufficiently  dense  sampling  grid  which  depends  on  the  known  support  constraint.  We  will 
not  consider  this  aspect  of  the  problem.  For  a  discussion  of  this  question,  see  [9],  In  the  case  where  the 
autocorrelation  function  is  known  directly,  then  the  support  of  jjm.nj  can  be  estimated. 
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Chapter  3 

Previous  Results  in  the  Phase 
Retrieval  Problem 

The  earliest  results  on  reconstruction  from  Fourier  transform  magnitude  for  sig¬ 
nals  with  given  support  relate  to  the  one-dimensional  case.  These  results  are  reviewed 
in  this  chapter,  both  because  they  provide  some  insight  into  the  general  problem  of 
phase  retrieval  and  also  because  solving  one-dimensional  problems  turns  out  to  be  an 
integral  part  of  several  proposed  algorithms  for  two-dimensional  phase  retrieval. 

The  second  part  of  this  chapter  discusses  the  two-dimensional  problem.  A  review  of 
the  known  results  on  the  uniqueness  and  characterization  of  the  solution  is  presented 
followed  by  a  review  and  critique  of  several  algorithms  which  have  been  proposed  for 
reconstruction  of  images  from  the  Fourier  transform  magnitude. 

3.1  Reconstruction  of  One-Dimensional  Signals  from 
Known  Support  and  Magnitude 

The  earliest  reference  to  the  one-dimensional  case  seems  to  be  a  paper  by 
Akutowicz  '11'  which  provided  a  complete  characterization  of  the  problem  of  one¬ 
dimensional  continuous  reconstruction  from  Fourier  transform  magnitude.  This  rela- 


tionship  was  rediscovered  simultaneously  by  Walther  [12]  and  Hofstetter  [13]. 

While  Akutowicz,  Walther  and  Hofstetter  considered  the  reconstruction  of  a  con¬ 
tinuous  one-dimensional  signal  from  its  Fourier  transform  magnitude,  more  recently 
Hayes  [9j  has  considered  the  discrete  one-dimensional  case.  The  situation  of  a  discrete 
one-dimensional  signal  can  be  derived  in  a  manner  similar  to  the  argument  followed 
by  Akutowicz,  Walther  and  Hofstetter.  His  development  is  given  here  while  stressing 
some  points  which  will  become  important  in  our  later  discussion. 

In  the  section  below,  we  discuss  virtually  simultaneously  the  problems  of  character¬ 
izing  the  number  of  solutions  to  the  one-dimensional  phase  retrieval  problem  and  an 
algorithm  for  producing  such  a  reconstruction.  This  is  in  contrast  with  our  discussion  of 
the  two-dimensional  phase  retrieval  problem  later  in  this  chapter,  where  reconstruction 
algorithms  are  developed  at  a  distance  from  uniqueness  considerations. 

The  object  of  the  algorithm  is  to  find  all  signals  of  extent  [0,  N\  which  have  the 
Fourier  transform  magnitude  \X(e*r)\.  From  (2.10)  we  know  that  we  can  immediately 
calculate  the  autocorrelation  function  of  any  such  signal  via, 

F-*{|X(e'-)|*}  =  r[n[  (3.1) 

Since  the  autocorrelation  function  of  an  N  -f  1  point  signal  has  support  [-W,  W]  the 
z-transform  of  r[n]  becomes 

*(*)=  £  rIn!2  "  (3-2) 

Tl=Z  —  X 

The  polynomial  associated  with  r[nj  is 

Vr{z)  =  tlr\N  ~  n]»"  (3-3) 

n=0 
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The  next  step  is  to  develop  the  product  decomposition  of  pr(z).  Since  r[nj  =  r*[  — n], 
pr(z )  is  its  own  mirror  polynomial.  Thus  if  pr(zo)  =  0,  then  p,.^-1)  =  0  also.  Further¬ 
more,  since  the  Fourier  transform  of  r[n]  is  non-negative,  the  zeros  z0  and  z q-1  will  be 
distinct,  or  occur  in  even  multiplicities.  Therefore,  p,(z)  can  be  expressed  as 

.v 

Pr(z)  =  All(z-Zi){l  -z\z)  (3.4) 

«=1 

where  A  is  positive.  Now  suppose  that  x[n]  with  associated  polynomial  p2{z)  is  one 
solution,  i.e.,  x\ n]  has  rinj  as  an  autocorrelation  function.  Recall  that 

Pr(z)  =  P*(Z)P*(Z)  (3-5) 

The  object  of  the  discussion  below  is  to  generate  a  p2(z )  that  satisfies  (3.5).  From 
the  set  of  zeros  in  (3.4),  one  zero  is  picked  from  each  pair  (z,-,2*-1)  and  a  polynomial 
with  such  zeros  is  generated  but  with  an  unknown  scale  constant  a , 

pAz)  =  aFI(2  -  *«•)  IK1  - zz))  (3-6) 

«€/ 

where  I  is  a  subset  of  ;1,  N].  The  mirror  polynomial  of  p,(z)  is  p2(z)  and  has  product 
decomposition  given  by 

P*(z)  =  <*'  IlC1  -  22>‘)  Il(2  -  *;)  (3-7) 

>ei  i<ii 

Now  the  product  of  p2{z)  and  p2{z)  will  equal  pr{z)  if  |a|2  =  A.  Thus,  a  can  be  picked 
to  be  any  complex  number  sucti  that  o  =  \Ta. 

From  the  procedure  above,  a  polynomial  associated  with  a  possible  solution  signal 
has  been  constructed.  Consider  two  different  polynomials  which  have  been  constructed 


using  the  algorithm  above.  They  will  have  some  zeros  in  common  and  some  zeros 


which  are  “flipped”  pairs.  Thus,  one  can  generate  alternate  polynomials  from  a  single 
construction  by  “zero- flipping”. 

At  this  moment  we  should  review  the  process  of  generating  the  associated  poly¬ 
nomial  of  a  solution  to  the  phase  retrieval  problem.  The  first  step  is  to  extract  the 
autocorrelation  function  from  the  Fourier  transform  magnitude  information.  From  this 
signal,  its  associated  polynomial  is  extracted.  The  zeros  of  this  polynomial  are  calcu¬ 
lated  and  these  zeros  occur  in  N  pairs.  From  each  of  these  N  pairs  one  member  is 
chosen  and  thus  determine  the  zeros  of  the  solution  associated  polynomial.  Since  one 
can  pick  either  member  from  each  pair,  a  total  of  2A  possible  zero  sets,  and  corre¬ 
sponding  polynomials  can  be  constructed.  The  final  step  is  to  find  an  appropriate  scale 
constant  for  the  solution  associated  polynomial.  This  scale  constant  has  a  prescribed 
magnitude  but  arbitrary  phase. 

If  the  desired  solution  must  be  real  then  zeros  chosen  for  a  solution  associated 
polynomial  must  include  complex  conjugate  pairs,  thereby  decreasing  the  number  of 
solutions  from  2V  possibly  down  to  2,V/'2  if  all  zeros  are  complex.  The  scale  constant 
is  also  restricted  to  being  real,  so  only  a  sign  ambiguity  instead  of  a  phase  ambiguity 
remains. 

Once  the  appropriate  associated  polynomial  has  been  specified  we  can  label  as  a 
solution  any  signal  which  is  associated  with  that  polynomial.  For  example,  if  the 
extracted  polynomial  is  given  by 

.V 

pAz)  =  Y.  a"z"  (3-g) 

n  =  0 
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then  a  possible  solution  is 


x[n]  =  a.v_n  for  0  <  n  <  N 
=  0  elsewhere 

In  fact  any  signal  of  the  form 


(3.9) 


x\n  +  k\=  a_v— n  for  0  <  n  <  N 

(3.10) 

=  0  elsewhere 
for  any  integer  A:  is  a  valid  solution. 

Although  the  algorithm  presented  in  the  previous  discussion  is  rather  simple  and  at 
first  sight  does  not  seem  to  deserve  the  painstaking  attention  which  has  been  directed 
to  it,  we  will  find  that  the  two-dimensional  problem  can  be  attacked  in  much  the  same 
way;  the  only  difference,  and  unfortunately  it  is  a  big  difference,  is  that  extracting 
the  appropriate  zeros  of  the  polynomial  associated  with  a  solution  will  be  much  more 
involved. 


3.2  Reconstruction  of  Two-Dimensional  Signals  from 
Known  Support  and  Magnitude 


3.2.1  Theoretical  Results 

Although  the  one-dimensional  reconstruction  from  Fourier  transform  magnitude 
given  a  finite  support  has  been  completely  solved,  the  two-dimensional  problem  has 
proven  to  be  both  more  challenging  and  more  interesting.  The  key  to  characterizing 
the  general  reconstruction  from  Fourier  transform  magnitude  problem  is  to  note  how  the 
multiplicity  of  solutions  in  the  one-dimensional  problem  is  introduced.  This  multiplicity 
depends  on  the  ability  to  decompose  the  polynomial  associated  with  a  solution  into  a 


product  of  lower  order  polynomials,  in  fact,  polynomials  of  degree  1.  It  is  this  property 
that  allows  for  the  “zero-flipping”  described  above,  resulting  in  several  solutions.  We 
will  discuss  the  fact  that  most  two-dimensional  signals  do  not  allow  such  zero-flipping; 
hence,  the  solution  to  the  phase  retrieval  problem  becomes  essentially  unique. 

The  original  insight  on  the  relationship  between  factorability  of  the  z-transform  and 
uniqueness  of  the  Fourier  transform  magnitude  reconstruction  was  presented  by  Bruck 
and  Sodin  [lj,  and  later  formalized  by  Hayes  [2]. 

Consider  a  known  Fourier  transform  magnitude  |X(e,u,e;”)|.  We  want  to  find  a 
signal  x[m,  nj  of  support  [0,A/]  x  [0,  N]  with  the  given  Fourier  transform  magnitude. 
From  the  Fourier  transform  magnitude  we  can  calculate  the  autocorrelation  function 
r[m,  nj  which  must  have  support  [~M,  M]  x  [-JV,  JVj.  The  polynomial  associated  with 
r[m,  n]  must  be  of  the  form, 

pr(w,z)  =  p,{w,z)p,{w,z)  (3.11) 

Thus,  pr(w,  z )  must  satisfy  two  conditions:  first,  it  must  be  a  polynomial  factorable  into 
two  smaller  polynomials  of  degree  ( M,N )  each.  Second,  each  of  the  factors  must  be 
mirror  polynomials  of  each  other.  Any  factorization  of  pr(w ,  z )  which  satisfies  these  two 
conditions  corresponds  to  a  valid  signal  z[m,  nj  which  satisfies  the  original  information. 

An  important  question  is  whether  there  is  only  one  factorization.  To  answer  this 
question  it  is  necessary  to  look  at  the  factorability  of  pt(w,z)  itself.  Following  Hayes 
'9'  define  the  following  equivalence  class, 

yjm.n)  ~  x[m,  nj  if  y[m,nj  =  e^xjfc,  ±m,k2±  n]  (3-12) 

for  some  integer  pair  {ki,k2)  and  some  6.  Thus,  y[m, n]  and  x[m, nj  are  equivalent  if 
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they  are  related  by  a  linear  shift,  a  phase  change,  and/or  a  rotation  by  180  degrees. 
From  the  properties  of  associated  polynomials,  we  see  that  x[m,  n]  ~  y[m,  nj  is  the 
same  as  stating  that  py(w,z)  =  exp(;0)p,(tu,  z)  or  py(w,z)  =  exp(/0)pJ(to,  z).  Using 
this  definition,  we  can  state  a  uniqueness  theorem  for  reconstruction  from  Fourier 
transform  magnitude, 

Theorem  3.1  If  x[m,  nj  has  an  irreducible  associated  polynomial,  then  all  other  y[m,  nj 
which  have  the  same  Fourier  transform  magnitude  must  be  equivalent  to  i[m,n). 

The  complete  proof  of  this  statement  is  given  in  [9]  for  the  case  of  real  z[m,  nj.  The 
case  for  complex  x[m,  n]  can  be  shown  in  a  similar  manner. 

It  has  been  shown  that  factorable  polynomials  form  a  zero-measure  set  in  the  space 
of  two-dimensional  signals  [14].  A  subsequent  contribution  by  Sanz  et  al.  [15]  showed 
that  not  only  do  signals  with  a  factorable  z-transform  form  a  zero-measure  set,  but  they 
also  form  a  non-dense  set  in  the  space  of  two-dimensional  signals.  From  these  results  we 
can  conclude  that  most  two-dimensional  signals  have  an  associated  polynomial  pt{w,  z ) 
which  is  irreducible.  Moreover,  the  polynomial  remains  irreducible  if  the  signal  is 
perturbed  by  a  small  amount. 

Several  comments  are  in  order  regarding  Theorem  3.1.  We  want  to  note  first  of 
all  that  a  finite  extent  constraint  must  be  placed  on  y\m,  nj,  i.e.,  in  the  Hayes  (and 
Bruck  and  Sodin)  proof,  it  is  implicitly  assumed  that  yjm,  n]  is  a  finite  extent  signal. 
Otherwise,  “infinite  order  polynomials”  do  not  satisfy  a  unique  factorization  property 
of  Theorem  2.1.  This  observation  implies  that  some  knowledge  of  the  support  of  the 
signal  is  needed  in  order  to  assure  a  well-behaved  uniqueness  result.  For  example,  if 


the  support  of  z(m,n]  is  [0, 3]2,  there  may  be  only  essentially  one  signal  of  finite  extent 
with  the  same  Fourier  transform  magnitude  as  x[m,  n]  but  there  may  be  many  signals 
of  extent  [0, 100]2  with  almost  the  same  Fourier  transform  magnitude  as  x[m,  n]. 

Under  some  conditions,  it  may  be  possible  to  specify  the  z[m,  n]  with  a  given  Fourier 
transform  magnitude  to  a  greater  extent  than  just  the  fact  that  it  is  a  member  of 
the  equivalence  class  defined  above.  For  example,  if  the  signal  is  composed  of  image 
intensities,  then  the  sign  ambiguity  disappears.  Also,  if  the  support  of  the  signal  is 
known  to  be  such  that  it  is  not  invariant  to  a  180  degree  rotation,  for  example  a 
triangular  support,  then  the  rotation  ambiguity  is  removed.  However,  if  the  support  is 
symmetric,  for  example,  square,  or  only  bounded  by  a  symmetric  region,  then  we  cannot 
distinguish  between  the  original  image  and  the  image  reversed.  This  ambiguity  is  all 
that  is  left  of  the  multiplicity  of  solutions  which  was  observed  in  the  one-dimensional 
case  that  was  due  to  the  “zero-flipping”;  either  no  zeros  are  flipped,  or  they  all  are. 
The  linear  shift  ambiguity  itself  is  seldom  of  importance. 

The  continuous  two-dimensional  case  has  only  been  studied  in  depth  recently.  An 
early  discussion  is  due  to  Huiser  and  Van  Toorn  [16]  who  first  studied  the  question. 
More  recently,  Sanz  et  al.  [17]  showed  that  the  non-uniqueness  of  the  case  studied 
by  Huiser  was  due  strictly  to  the  fact  that  the  Laplace  transform  Huiser  considers  is 
factorable,  and  in  fact  showed  that  multidimensional  continuous  signals  with  Laplace 
transforms  which  are  non-factorable  entire  functions  can  be  uniquely  reconstructed 
from  the  magnitude  of  the  Fourier  transform. 


3.2.2  Study  of  Previous  Algorithms  for  Phase  Retrieval 


Fienup  Algorithm 

The  development  of  algorithms  for  reconstruction  of  two-dimensional  signals  from 
Fourier  transform  magnitude  preceded  the  establishment  of  conditions  for  uniqueness. 
The  earliest  family  of  algorithms  were  of  an  iterative  nature,  pioneered  by  Fienup 
[18,19’.  These  algorithms  are  based  on  work  by  Gerchberg  and  Saxton  [20],  and  since 
the  introduction  by  Fienup  have  been  studied  by  Hayes  [2],  Sanz  and  Huang  [21],  Levi 
and  Stark  [22],  and  Won  et  al.  [23],  among  others.  In  this  section  we  describe  and  study 
three  versions  introduced  by  Fienup.  The  simplest  algorithm  is  the  Gerchberg-Saxton- 
Fienup  (GSF)  algorithm  which  iteratively  imposes  the  known  support  in  the  spatial 
domain  and  the  known  Fourier  transform  magnitude  in  the  frequency  domain.  The 
second  algorithm  considered  is  the  Gerchberg-Saxton-Fienup  algorithm  with  positivity 
constraint  (GSF-P)  which  forces  each  estimate  to  be  non-negative  as  well  as  of  the  given 
support.  Although  positivity  is  not  a  requirement  for  uniqueness  of  reconstruction,  it  is 
hoped  that  using  such  additional  information  will  aid  in  speed  of  convergence.  The  final 
algorithm  considered  is  the  hybrid  input-output  (HIO)  algorithm  developed  by  Fienup 
which  is  described  below.  The  object  of  this  section  is  to  describe  these  algorithms  and 
evaluate  their  performance. 

The  algorithms  begin  with  an  initial  estimate  which  can  be  derived  either  from  a 
prion  knowledge  of  the  signal,  i.e.,  a  low  resolution  image,  or  in  the  case  of  no  other  a 
priori  knowledge,  from  an  image  with  either  a  random  phase  component  which  satisfies 
the  Fourier  transform  magnitude  constraint,  or  from  an  image  with  random  coefficients 


which  satisfies  the  support  constraint.  In  our  discussion,  we  take  the  last  approach. 


Given  an  estimate  x,[m,  n]  of  x[m,  n]  which  satisfies  the  support  constraints,  an 
auxiliary  signal  x,[m,  n]  is  calculated  by  the  three  algorithms,  where  x,[m,  n]  has  the 
correct  Fourier  transform  magnitude  and  the  same  phase  as  x,[m,  nj, 

x,[m,n]  =  F “ 1  { | X( p ) | eJ 0,9 * •v‘ ^  }  V  (m,n)  (3.13) 

where  Xt{e’u,e>v)  is  the  Fourier  transform  of  x,[m,  n).  However,  in  general,  x,[m,  n]  will 
not  satisfy  the  known  support  and/or  positivity  constraints.  The  difference  between 
the  three  algorithms  is  how  the  next  estimate  x.+Jm,  nj  is  calculated  from  x,[m,  n]  and 
x.  jm,  nj. 

For  GSF,  x1+l[m,  n]  is  calculated  as, 

x,[m,  n]  (m,  n)  G  S 

x,+1(m,  nj  =  •  (3.14) 

0  elsewhere 

where  the  S  denotes  all  index  pairs  where  x[m,  n]  is  allowed  to  be  non-zero.  In  the 
GSF-P  algorithm,  xl+l(m,nj  is  derived  from, 

x,[m,  nj  ( m,n)  €  S  and  x,[m,  nj  >  0 
xt+1(m,  nj  =  '  (3.15) 

0  elsewhere 

Finally,  HIO  uses  both  the  previous  estimate  x,[m,nj  and  x,[m,n], 

x,!m,n]  (m,n)e$  and  x,[m,nj>0 

x.+i  [m.  n!  =  (3.16) 

x,jm,  nj  —  x,[m,  nj  elsewhere 

From  this  estimate,  x,+i[m,nj,  the  iteration  is  repeated.  We  note  that  implementa¬ 
tion  requirements  dictate  that  the  algorithm  use  a  sampled  Fourier  transform.  Moreover 
in  order  to  be  able  tc  impose  a  support  constraint,  we  need  to  ov^rsampie  the  Fourier 


transform  magnitude.  In  our  examples,  if  the  support  is  enclosed  in  an  N  by  N  “box” 
we  use  a  2 N  by  2 N  discrete  Fourier  transform  of  the  appropriately  zero-padded  image. 

There  is  considerable  disagreement  between  researchers  in  the  field  considering  the 
conditions  under  which  these  iterative  algorithms  converge  to  the  correct  reconstruc¬ 
tion.  It  can  be  shown  [19]  that  the  GSF  algorithm  decreases  at  each  iteration  the  mean 
squared  error  between  the  estimate  and  true  Fourier  transform  magnitude.  Thus,  the 
algorithm  must  converge  to  a  fixed  point  of  the  iteration.  This  convergence  does  not 
imply,  however,  that  the  error  will  decrease  to  zero. 

Fienup  [18,19]  has  reported  on  the  wide  applicability  of  the  iterative  algorithms, 
especially  the  HIO  algorithm.  On  the  other  hand,  Hayes  [2]  has  found  that  the  GSF 
algorithm  was  not  successful  in  reconstructing  several  original  images.  In  a  response, 
Fienup  [24]  has  noted  that  the  positivity  constraint,  which  is  used  in  the  GSF-P  and 
HIO  algorithms,  is  important  in  ensuring  convergence  to  the  original  image.  However, 
Sanz  et.  al  [21]  has  produced  an  example  where  the  GSF-P  algorithm  did  not  converge 
to  the  original  image  even  though  the  algorithm  makes  use  of  the  positivity  constraint. 

In  Appendix  A  we  present  a  study  of  the  convergence  properties  of  the  algorithms 
described  above.  A  Monte  Carlo  study  was  performed  on  a  large  number  of  randomly 
generated  4  by  4  and  8  by  8  images  to  assess  the  convergence  rate  of  each  algorithm 
presented.  The  results  of  this  study  are  summarized  in  Tables  3.1  and  3.2.  These  tables 
present  the  percent  of  trials  which  converged  to  a  solution  which  had  a  normalized  mean 


squared  error  (defined  in  Appendix  A)  of  either  less  than  10  2  or  less  than  10  4. 


Success  rate  of  convergence  (percent)  | 

support  type 

rectangular 

triangular 

cutoff 

10”2  10'4 

lO'2  10'4 

GSF 

21  9 

59  53 

GSF-P 

22  12 

60  55 

HIO 

22  12 

60  52 

Table  3.1:  Summary  of  algorithm  convergence  as  a  function  of  support  and  criterion 
for  4  by  4  image. 


Success  rate  of  convergence  (percent) 

support  type 

rectangular 

triangular 

cutoff 

10'2  10“4 

10-2  10'4 

GSF 

0  0 

44  40 

GSF-P 

0  0 

54  50 

HIO 

0  0 

54  48 

Table  3.2:  Summary  of  algorithm  convergence  as  a  function  of  support  and  criterion 
for  8  by  8  image. 

From  Tables  3.1  and  3.2  we  see  that  the  iterative  algorithms  perform  substantially 
different  for  the  case  of  symmetric  and  non-symmetric  supports.  For  the  case  of  a 
non-symmetric  support,  the  success  rate  was  on  the  order  of  40  to  60  percent  while  for 
symmetric  supports,  the  algorithms  succeeded  in  at  most  20  percent  of  the  trials.  The 
reason  for  this  difference  may  be  that  for  non-symmetric  supports,  there  is  no  rotation 
ambiguity  in  the  reconstruction,  i.e.,  one  cannot  rotate  the  image  by  180  degrees  and 
be  able  to  “fit”  the  image  in  the  same  support.  Moreover,  while  the  convergence  rate 
for  the  non-symmetric  support  remained  about  the  same  as  the  image  size  increased, 
the  iterative  algorithms  degraded  dramatically  for  8  by  8  square  images  compared  to 
4  by  4  square  images. 

We  must  make  the  final  observation  that  although  there  was  a  difference  in  the  con¬ 
vergence  rate  among  the  three  algorithms,  the  success  rate  among  the  three  algorithms, 
as  defined  in  Tables  3.1  and  3.2  is  about  the  same;  that  is,  the  frequency  with  which 
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the  application  of  any  of  the  algorithms  considered  resulted  in  a  good  reconstruction 
of  the  original  image  was  about  the  same. 

There  are  several  conclusions  we  can  derive  from  the  study  presented  above.  For 
the  family  of  symmetric  signals  considered,  we  can  surmise  that  the  algorithms  consid¬ 
ered  converged  to  the  correct  reconstruction  in  at  most  20  percent,  depending  on  the 
reconstructed  image  fidelity  cutoff  used  and  the  image  size.  As  the  image  size  increased 
the  success  rate  decreased.  In  the  case  where  a  nonsymmetric  region  of  support  was 
known,  the  iterative  algorithm  performed  substantially  better,  achieving  a  final  esti¬ 
mate  close  to  the  true  signal  in  about  50  to  60  percent  of  the  trials,  the  success  rate 
staying  constant  as  the  image  size  was  increased.  The  use  of  positivity  information  did 
not  aid  the  convergence  rate  substantially. 


Bates  Algorithm 


Bates  [25,26,27]  has  developed  an  algorithm  for  the  two-dimensional  reconstruction 
which  tries  to  estimate  the  phase  of  the  Fourier  transform  directly  from  samples  of  the 
Fourier  transform  magnitude.  For  simplicity,  consider  a  two-dimensional  signal  x[m,  n] 
which  has  a  region  of  support  L0,  X]‘.  Consider  the  samples  of  the  Fourier  transform 
at  frequencies  u  =  2nk/(N  +  1),  v  =  2nl/{N  +  1), 

XkJ  =  (3.17) 


Since  X{eJU,eJV)  is  the  Fourier  transform  of  an  [0,  signal,  it  must  satisfy  the  fol¬ 
lowing  “interpolation”  formula, 


X{e,'\ejr)  -  e;'NHi  33  31  sinc(- -  u)  smc{  ~~  -  v) 
*=o  i=o  N  +  1  ‘V  *  1 


(3.18) 


where 
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The  algorithm  begins  by  setting  the  phase  of  X0,o  to  zero.  As  Bates  points  out, 
this  is  not  a  restriction  to  the  algorithm  since  all  signals  which  differ  by  a  constant 
phase  have  the  same  Fourier  transform  magnitude.  The  algorithm  then  proceeds  to 
calculate  the  two  possible  approximations  to  the  phase  of  Xi.o  via  (3.22).  One  of  these 
phase  approximations  is  picked  arbitrarily.  Again  this  is  not  really  a  restriction  since 
the  only  effect  is  to  possibiy  rotate  the  resulting  reconstruction  by  180  degrees.  The 
two  possible  values  for  the  phase  of  X0A  are  then  calculated  via  (3.23).  At  this  point 
both  estimates  for  the  phase  of  X0.i  must  be  kept  since  all  the  degrees  of  freedom  have 
been  used  in  setting  the  phase  of  X0.0  and  The  algorithm  then  proceeds  by  using 
both  possible  values  for  the  phase  of  Xo.i  and  (3.22)  to  calculate  four  possible  values 
for  the  phase  of  X^. 

At  this  point  the  phase  of  X0,o  and  X0.i  is  known  approximately,  there  are  two 
possible  phase  values  for  Xi,0  and  four  possible  phase  values  for  X^,.  The  next  step  is 
to  remove  the  ambiguity  in  the  phase  of  Xlml  by  re-calculating  the  phase  of  XlA  in  an 
independent  manner,  namely  via  the  knowledge  of  the  phase  of  X1-0.  Since  the  phase 
of  X^n  has  been  uniquely  specified,  we  can  use  (3.22)  to  calculate  only  two  possible 
values  of  the  phase  of  XlA.  Bates  maintains  that  the  two  sets  of  possible  values  for  the 
phase  of  X^i  derived  via  the  two  different  paths 

X„.o  — *  X».0  — >  (3-24) 

and 

X0.o  — *  X0.i  —*  XlA  (3.25) 

will,  in  general,  have  oniy  one  value  in  common,  or  more  precisely  only  one  pair  which 


almost  matches.  Using  this  value,  the  phase  of  is  thus  uniquely  specified.  Tracing 
our  steps  back,  we  can  then  get  a  unique  phase  estimate  for  .Jfo.i-  The  algorithm  then 
proceeds  to  calculate  other  values  of  phase  by  working  sequentially  away  from  X0.o- 

Obviously,  this  algorithm  makes  many  assumptions  which  make  it  a  questionable 
solution  to  the  phase  retrieval  problem  although  the  authors  have  had  some  success  in 
implementation  [26,27].  First  of  all,  the  whole  algorithm  is  predicated  on  the  linear 
assumption  made  in  (3.20);  this  is  rarely  valid.  Second,  as  the  authors  suggest,  this 
assumption  will  cause  the  two  sets  of  phase  values  for  Xla  not  to  have  any  points  in 
common,  necessitating  the  use  of  a  “closest”  criterion.  As  soon  as  one  error  is  made 
in  choosing  a  phase  value,  all  later  phase  estimates  will  be  incorrect.  Even  if  no  error 
is  made  in  choosing  the  “closest”  match,  higher  frequency  phase  estimates  will  be 
increasingly  poor. 

Because  the  phase  corresponding  to  higher  frequency  terms  get  progressively  poorer, 
we  expect  that  the  algorithm  would  be  able  to  extract  with  some  accuracy  the  lower 
frequency  content  of  the  image.  In  fact,  this  is  the  observed  phenomenon;  the  recon¬ 
struction  seems  to  be  a  severely  blurred  version  of  the  original  [23].  More  recently,  some 
effort  has  been  made  in  using  the  Bates  algorithm  to  generate  an  initial  estimate  for 
the  Fienup  algorithm.  For  the  case  where  the  image  has  mostly  low  frequency  content, 
the  Bates  algorithm  gives  an  adequate  initial  estimate  for  the  Fienup  algorithm  [28]. 

Canterakis  Algorithm 

Recently,  a  third  algorithm  has  been  proposed  by  Canterakis  [29].  This  algorithm 
differs  from  the  algorithms  presented  so  far  in  that  it  is  guaranteed  to  yield  an  exact 
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solution  to  the  phase  retrieval  problem.  The  algorithm  can  be  described  as  a  conver¬ 
sion  of  the  two-dimensional  reconstruction  problem  to  a  one-dimensional  reconstruction 
problem.  Thus,  it  is  important  to  keep  in  mind  the  discussion  of  Section  3.1  in  evalu¬ 
ating  this  algorithm.  For  simplicity,  suppose  that  the  support  of  the  unknown  signal 
z'm,  nj  is  [0,  JVj2 .  Consider  now  the  one-dimensional  signal  z[n]  given  by 

xm  -I-  (2N  -f  l)n]  =  x{m,  nj  for  0  <  m,n  <  N  (3.26) 

=  0  elsewhere  (3.27) 

This  can  be  viewed  as  attaching  N  zeros  to  each  row  of  x[m,  nj  and  concatenating  the 
rows.  It  is  easily  shown  that  f  [n] ,  the  autocorrelation  function  of  x[n],  satisfies 

f[m  +  (2iV  +  l)n]  =  r[m,  n]  (3.28) 

where  rim,  n]  is  the  autocorrelation  function  of  x[m,  nj.  Note  moreover  that  (3.28)  is  a 
reversible  transformation;  f  [nj  is  generated  from  r[m,  n]  by  concatenating  the  rows  of 
rjm,  nj.  The  Canterakis  algorithm  consists  of  first  calculating  r[m,nj  from  the  known 
Fourier  transform  magnitude  via  the  relationship  of  (2.13).  Then  fjnj  is  formed  and 
all  solutions  to  the  one-dimensional  reconstruction  problem  are  calculated.  Since  z[nl 
is  an  2N2  4-  2 N  -f  1  point  signai,  there  are  about  2V’  solutions  to  the  one-dimensional 
problem  posed,  taking  into  account  the  fact  that  zV  is  assumed  to  be  real.  Each 
candidate  i  n  is  tlien  checked  to  see  if  it  can  correspond  to  a  x:m.  ni;  nameiy,  whether  it 
consists  of  the  alternation  of  non-zero  and  zero  samples  consistent  with  a  “transformed” 
two-dimensional  signal.  For  each  successful  x[nj,  the  corresponding  z[m,  nj  is  formed, 
resulting  in  a  solution  to  the  two-dimensional  reconstruction  problem. 


It  is  clear,  as  the  author  pointed  out,  that  this  algorithm  requires  exponentially  in¬ 


creasing  computations  as  the  signal  size  increases.  Because  there  are  on  the  order  of  2-v* 
solutions  to  the  associated  one-dimensional  problem,  the  computational  load  increases 
rapidly  as  N  increases.  Clearly  this  algorithm  cannot  begin  to  be  applied  to  signals  of 
reasonable  size.  We  should  note,  however,  that  the  Canterakis  algorithm  does  provide 
an  interesting  perspective  on  reconstruction  from  magnitude.  One  point  is  that  the 
phase  retrieval  problem  is  posed  as  that  of  extracting  x [m,  n]  from  its  autocorrelation 
function  instead  of  extracting  x[m,  nj  from  the  Fourier  transform  magnitude  directly. 
The  second  interesting  point  is  that  the  Canterakis  method,  as  we  will  see  later,  is 
an  implicit  attempt  at  factoring  the  polynomial  associated  with  the  autocorrelation 
function  of  x[m,  n].  This  is,  in  fact,  the  specific  approach  we  will  pursue  in  Chapters  4 
and  5. 

Deighton,  Scivier  and  Fiddy  Algorithm 

Recently  Deighton,  Scivier  and  Fiddy  [30j  have  developed  an  algorithm  reminis¬ 
cent  of  the  Bates  algorithm.  This  algorithm  however  does  not  make  any  approxima¬ 
tions,  and  therefore,  like  the  Canterakis  algorithm,  yields  an  exact  reconstruction. 

Suppose  that,  as  in  the  discussion  of  the  Bates  algorithm,  one  seeks  to  calculate  the 
phase  difference  between  the  Fourier  transform  samples  Xtj  and  Xk+ij.  Let  ij[n]  be 
defined  by 

.v 

x,\n)  =  i[m,  n]e_;^+?  (3.29) 

m— 0 


where  again,  x\m,  nj  has  support  [0,  fVj2,  and  let  the  Fourier  transform  of  i/jnj  be 


denoted  by  X/(e;<).  Using  x/[nj  above,  we  can  express  X*,/  and  Xi+1,j  as 


XkJ  =  *,(«*£) 


Xk+u  =  X,(e  -v-m  ) 


(3.30) 

(3.31) 


Thus  the  Fourier  transform  of  x/[n]  is  the  strip  in  the  (u,  v)  plane  along  v  = 
(2 ?r/)/(JV  4-  1).  Since  jX(e;“,  e;,')|  is  known  for  all  (u,t»),  the  Fourier  transform  mag¬ 
nitude  of  X|[nj  is  also  known  for  all  frequencies.  The  signal  xf[n]  has  support  [0,  N\ 
and  the  results  on  solving  the  one-dimensional  phase  retrieval  problem  can  be  used  to 
generate  a  set  of  at  most  2V  signals  of  which  one  must  be  X/[n].  By  calculating  the 
Fourier  transform  phase  of  each  of  these  possible  signals,  the  possible  phase  difference 
between  X*./  and  Xk+u  is  thus  restricted  to  be  one  of  at  most  2-v  different  possibilities. 
From  the  same  process,  2V  different  possible  values  for  the  phase  difference  between 
Xk.i+i  and  Xk-ri.i+i  can  be  generated.  Combining  the  two  results  above,  we  have  22V 
possible  phase  differences  between  Xk.i  and  Xk+u+i-  By  using  the  Xkj  — ►  Xjt.i+i  — »  Xk.i 
path,  we  get  another  2:A  possible  phase  differences.  Comparing  the  two  lists,  we  expect 
that  only  two  phases  will  be  in  both  lists,  the  true  phase  and  the  phase  of  the  reversed 
image.  Although  no  proof  has  been  presented,  experimentally  it  has  been  found  that 
only  those  two  phase  values  will  be  in  common  in  both  lists. 

The  Deighton  algorithm  proceeds  as  follows:  first,  the  phase  of  X0.0  is  arbitrarily  set 
to  zero.  Via  the  phase  matching  process  just  described  the  phase  of  Xu.!,  X10  and  Xt.! 
are  calculated.  Moreover,  the  phase  of  X, ».*,  Xi .*  and  Xk.\  can  be  calculated  for  all  k. 
Suppose  that  the  phase  of  Xt.;  is  now  desired  for  k  >  2.  Using  the  same  procedure,  all 
2  V  solutions  for  x;  n  are  calculated.  The  one  solution  which  agrees  with  the  previously 


as 


Figure  3.1:  Phase  differences  to  be  used  in  Deighton  algorithm. 


calculated  phase  difference  between  X0.2  and  AT^;  is  then  used  to  calculate  the  phase 
of  Xk.2  for  all  other  values  of  k.  This  way,  the  rest  of  the  phase  values  of  Xk.i  can  be 
calculated. 

The  most  computationally  demanding  part  of  the  Deighton  algorithm  is  the  ini¬ 
tial  matching  of  phases.  Figure  3.1  displays  the  number  of  possible  phase  differences 
among  the  four  Fourier  components  Xo.o.  X0A,  XL0  and  XlA.  The  characterization 
of  the  number  of  phase  calculations  needed  turns  out  to  involve  a  tradeoff  between 
computation  and  storage.  Consider  first  trying  to  find  a  match  for  the  phase  difference 
between  X,.0  and  .VL1  in  the  manner  just  discussed.  Then,  the  algorithm  needs  to 


compare  two  lists  with  2 


.v-i-'.w: 


elements  each.  Thus  the  lists  grow  at  a  rate  of  about 


.  Such  a  list  becomes  too  large  to  be  accommodated  in  main  memory  (6  Mbytes  of 


available  memory)  for  images  larger  than  12  by  12.  Since  we  will  be  dealing  with  larger 
images,  it  is  imperative  that  the  memory  requirements  of  the  Deighton  algorithm  be 
reduced  for  comparison  purposes.  This  can  be  done  by  finding  a  match  for  the  phase 
difference  between  Xi.0  and  X0,i  instead  of  X0.o  and  Xj.i.  This  allows  us  to  use  a  list 
of  size  approximately  221-V/2J  so  we  could  consider  images  of  size  up  to  about  18  x  18. 
This  storage  reduction  comes  at  a  price  however,  since  now  the  number  of  comparisons 
required  is  about  22-v  instead  of  23-V/'2. 

To  the  author’s  knowledge,  of  all  the  algorithms  so  far  published  for  phase  retrieval 
which  is  guaranteed  to  yield  the  correct  solution,  the  Deighton  algorithm  is  the  most 
computationally  efficient.  In  order  to  be  able  to  compare  the  performance  of  our 
algorithm  with  the  Deighton  algorithm,  the  Deighton  algorithm  was  implemented  using 
the  second  method  described  here.  In  Figure  3.2  we  show  the  number  of  CPU  seconds 
on  a  Vax-750  required  to  reconstruct  an  image  with  support  [0 ,7V  -  l]2  as  a  function 
of  N. 

The  severe  storage  or  computational  requirements  of  Deighton  algorithm,  especially 
its  quick  growth  with  image  size,  is  directly  related  to  the  combinational  approach  used 
of  trying  all  possible  solutions.  Since  the  resources  required  grow  exponentially,  one 


quickly  runs  into  a  storage  or  computational  barrier. 


Chapter  4 


Phase  Retrieval  via  Bivariate 
Polynomial  Factorization 


We  have  so  far  discussed  several  approaches  to  reconstruction  from  Fourier  trans¬ 
form  magnitude.  The  algorithms  of  Fienup  and  Bates  are  computationally  efficient  and 
easily  applied.  However,  there  is  no  assurance  that  using  these  algorithms  will  yield  a 
solution  to  the  phase  retrieval  problem.  We  have  also  considered  two  algorithms  which 
provide  a  closed  form  solution  to  the  reconstruction  problem.  The  first  algorithm  stud¬ 
ied.  developed  by  Canterakis,  rapidly  becomes  prohibitively  expensive  as  the  image 
size  is  increased,  or  equivalently,  as  the  image  resolution  is  increased.  The  second 
closed  form  algorithm,  studied  by  Deighton  et  al.,  is  more  attractive  computationally. 
However,  the  required  computational  and  storage  load  also  increase  very  rapidly  as  a 
function  of  image  resolution. 

As  described  in  the  last  chapter,  the  results  of  Bruck  and  Sodin,  Hayes,  and  Sanz 
point  to  the  intimate  relationship  between  conditions  for  uniqueness  of  reconstruction 
from  Fourier  transform  magnitude  and  the  factorability  of  polynomials.  Moreover,  we 
found  that  a  general  procedure  for  solving  the  one-dimensional  phase  retrieval  problem 
consists  of  factoring  the  polynomial  associated  with  the  autocorrelation  function.  The 


approach  we  introduce  in  this  chapter  in  order  to  generate  a  new  algorithm  for  phase 
retrieval  is  to  view  the  relationship  between  phase  retrieval  and  polynomial  factorization 
as  also  pointing  to  a  method  for  solving  the  two-dimensional  phase  retrieval  problem. 
We  will  also  find  that  this  approach,  although  explicitly  expressed  here,  is  the  implicit 
basis  for  the  Canterakis  algorithm  described  earlier. 

If  the  coefficients  of  the  unknown  signal  are  known  to  be  rational  numbers,  then  the 
theory  of  integer  polynomial  factorization  becomes  applicable.  This  is  the  approach 
considered  by  Berenyi,  Deighton,  and  Fiddy  [31].  There  are  several  difficulties  with 
the  rational  coefficient  assumption.  First  of  all,  an  infinitesimal  amount  of  noise  will 
render  the  algorithm  useless.  Second  of  all,  algorithms  for  factoring  polynomials  are 
very  computationally  intensive.  In  Berenyi’s  paper,  the  largest  image  considered  was 
a  6  by  6  image. 

One  novel  idea  of  this  thesis  is  to  consider  algorithms  for  factoring  polynomials 
over  the  reals  or  complex  numbers.  This  makes  the  problem  more  robust  to  noise  and 
also  applicable  to  the  case  where  the  signal  to  be  reconstructed  has  complex  values. 
The  final  result,  described  in  the  next  chapter,  is  a  new  factorization  algorithm  which 
leads  to  a  closed  form  solution  to  the  phase  retrieval  problem  which  is  much  more 
computationally  efficient  than  the  Canterakis  or  Deighton  algorithms.  As  a  result,  the 
reconstruction  of  images  of  sizes  up  to  25  by  25  becomes  feasible. 

In  this  chapter  we  first  pose  the  two-dimensional  phase  retrieval  problem  in  anal¬ 
ogy  to  our  previous  one-dimensional  presentation  of  Section  3.1  and  emphasize  how 
factorization  of  the  resulting  bivariate  polynomial  leads  to  a  solution  of  the  phase  re¬ 
trieval  problem.  This  leads  to  the  topic  of  this  chapter,  which  is  the  discussion  of 


several  known  algorithms  for  factoring  polynomials  in  two  variables  over  the  complex 
numbers.  In  order  to  do  this  we  first  introduce  some  results  from  algebraic  function 
theory.  We  follow  this  by  a  discussion  of  the  algorithms  of  Kronecker  and  Kaltofen  for 
factoring  bivariate  polynomials.  In  the  next  chapter  we  introduce  the  main  result  of 
this  thesis,  which  is  a  new  algorithm  for  factoring  polynomials  in  two  variables  over 
the  reals  or  complex  numbers  and  its  application  to  the  phase  retrieval  problem. 

4.1  Phase  Retrieval  as  Polynomial  Factorization 

We  have  already  found  that  the  polynomial  associated  with  the  autocorrelation 
function  is  always  factorable  into  the  product  of  px(w,z),  the  polynomial  associated 
with  the  unknown  image  x{m,n ]  and  the  mirror  polynomial  px(w,z), 

p,{w,z)  =  p,(w,z)p,(w,z)  (4.1) 

Moreover,  if  px(vu,z)  is  not  factorable,  then  px(w,z)  is  not  factorable  either,  and  thus 
(4.1)  is  the  only  uon-trivial  product  decomposition  of  pr(w,z).  Thus,  if  px(w,z )  is  not 
factorable,  we  can  generate  x'm,  n j  from  rim,  n]  by  factoring  pr{w,  z)  into  its  irreducible 
components.  This  observation  is  the  basis  for  the  result  by  Bruck  and  Sodin,  and  Hayes. 
From  (4.1)  we  see  the  possible  utility  of  algorithms  for  factoring  bivariate  polynomials; 
namely,  if  we  can  factor  pr[w,z)  then  we  will  retrieve  p,{w,z)  or  px(w,z),  and  from 
p,{w,  z )  or  pAw,  z)  we  can  extract  x'm,  nj  or  x!-m,  -nj. 


4.2  Zeros  of  Bivariate  Polynomials 


4.2.1  Motivation 

In  the  search  for  factorization  mechanisms  for  polynomials  over  the  complex  num¬ 
bers,  zeros  of  such  polynomials  play  a  crucial  role.  The  importance  of  zeros  is  essentially 
due  to  the  property  that  the  set  of  zeros  of  a  polynomial  is  equal  to  the  union  of  the 
set  of  zeros  of  its  factors.  For  example,  the  set  of  zeros  z,  of  a  univariate  polynomial 
is  equal  to  the  union  of  the  zeros  of  each  of  its  irreducible  factors,  namely  each  factor 
(z  —  z,).  The  same  situation  holds  for  polynomials  in  two  or  more  variables.  Unlike 
univariate  polynomials,  however,  the  zeros  of  polynomials  in  two  variables  cannot  be 
broken  up  into  independent  pieces  consisting  of  an  isolated  zero  each.  If  the  polynomial 
is  reducible,  its  set  of  zeros  can  be  broken  up  into  a  union  of  zeros  corresponding  to  its 
irreducible  factors.  The  main  difficulty  is  that  performing  such  a  dissection  is  a  much 
more  intricate  matter  in  the  bivariate  case  than  the  univariate  case.  For  this  reason 
the  discussion  below  is  concerned  with  developing  some  of  the  properties  of  zeros  of 
polynomials  in  two  variables. 

4.2.2  Algebraic  Functions 

We  already  discussed  that  there  are  complex  numbers  z,  associated  with  a  polyno¬ 
mial  p(z),  called  the  zeros  of  p(z),  such  that  p(z,)  =  0.  We  also  noted  that  the  number 
of  such  zeros  is  finite.  In  fact,  there  are  at  most  deg(p)  zeros  of  p(z). 

Bivariate  polynomials  have  a  much  more  interesting  set  of  zeros.  In  fact  a  large 
branch  of  mathematics  is  concerned  with  the  study  of  the  zeros  of  polynomials  in  two 


variables.  A  very  readable  discussion  of  this  topic  is  contained  in  [10].  Most  of  the 

development  below  is  based  on  this  reference.  In  order  not  to  stray  too  far  from  our 

main  subject,  some  of  the  relevant  results  are  rederived  here  in  a  restricted  yet  adequate 

form,  but  without  the  necessity  of  introducing  certain  peripheral  subjects. 

Recall  that  a  zero  of  the  polynomial  of  degree  (Af,  N) 

\t  .v 

p(w>z)  =  Y,  zL  am.nWmZn  (4.2) 

m— 0  n=0 

is  a  pair  of  complex  numbers  (w,,  z,)  such  that  p(u/,,  z,)  =  0.  However,  unlike  the  one¬ 
dimensional  case,  the  number  of  such  zeros  is  not  finite.  In  fact,  as  shown  below,  for 
any  complex  number  z  we  can  always  find  a  w  such  that  (10,2)  is  a  zero  of  p(tu,  2). 

Consider  p(w,  z )  written  as  a  polynomial  in  w  with  coefficients  which  are  polyno¬ 
mials  in  2, 

p{w,z)  =  p0(2)  +  Pi(z)w  +  ■  •  •  +  pM(z)uru  (4.3) 

where 

Pji2)  =  Y,P,.kZk  (4.4) 

k- 0 

For  each  value  of  2,  the  equation  p{w,z )  =  0  defines  a  function  u>(z)  implicitly  by 

p(w{z),z)  =  0  (4.5) 

This  function  is  called  an  algebraic  function  because  for  all  2,  w(z)  satisfies  the  algebraic 
or  polynomial  equation  (4.5).  This  function  ui(z)  will  be  multi-valued;  specifically,  it 
will  normally  consist  of  the  M  distinct  finite  roots  of  (4.3) 

wl{z),w2{z),  ■  ■  ,wsl{z)  (4.6) 


These  are  called  the  branches  of  the  algebraic  function. 


It  is  straightforward  to  generate  these  roots;  namely,  each  p;(z)  in  (4.3)  is  evaluated 
at  z,  reducing  (4.3)  to  a  polynomial  in  xv  with  numerical  coefficients.  The  roots  of  this 
one-dimensional  polynomial  is  exactly  what  we  mean  by  (4.6). 

A  point  z  =  a  in  the  complex  z-plane  where  there  are  M  distinct  finite  zeros  of 
p(w,  a),  i.e.,  where  all  the  branches  of  u>(z)  are  finite  and  do  not  intersect,  is  called 
an  ordinary  point  of  the  algebraic  function  w(z).  A  point  where  one  of  the  branches 
becomes  infinite  or  two  branches  coalesce  into  a  single  root  is  called  a  singular  point. 
Therefore,  each  point  in  the  2-plane  is  either  an  ordinary  or  singular  point  of  w(z). 
The  case  where  two  branches  coalesce  into  a  single  root  is  also  called  a  branch  point. 
Thus,  a  branch  point  is  a  special  case  of  a  singular  point.  Singular  points  which  occur 
for  finite  values  of  z  will  be  called  finite  singular  points. 

Example:  Consider  the  polynomial, 

p(w.z )  =  l-z  +  z2-f(-l  +  z-  z2)w  +  (1  -  z)w2  (4.7) 

A  plot  of  w(z)  as  a  function  of  z  between  .5  and  1.5  is  shown  in  Figure  4.1.  In  this  case  there 
are  two  branches  since  the  degree  of  p(w.z)  in  w  is  2.  Thus  for  almost  all  z.  there  are  two 
values  of  w  such  that  p(w.z)  =  0.  One  singular  point  of  w(z )  occurs  at  z  =  1.  where  it  can 
be  seen  from  the  figure  that  one  branch  becomes  infinite  at  this  value  of  z.  A  branch  point  of 
w(z)  occurs  at  approximately  z  =  7.913.  w  =  2.01G47  where  the  two  branches  join.  For  real  z 
less  than  the  location  of  the  branch  point,  the  branches  of  w(z)  are  complex  conjugate  pairs. 
Figure  4.1  in  this  case  shows,  in  dotted  lines,  the  real  part  of  w  plus  or  minus  the  imaginary 
part  of  w. 


The  set  of  singular  points  can  be  characterized  by  the  following  two  theorems, 


Theorem  4.1  Let  p(w,  z)  be  a  polynomial  of  degree  (M,  N),  expressed  as  in  (4-8),  then 
all  w,(z )  evaluated  at  za  must  be  finite  unless  p.\/(zn)  =  0,  where  pu(z)  is  defined  in 
(4-3).  The  number  of  such  z0  is  at  most  N . 


Proof:  We  first  need  to  invoke  a  result  by  Cauchy  [32,  p.  123] , 


Lemma  4.1  All  the  zeros  of 

f(z)  =  a0  +  ayz  +  a2z2  +  ...  +  a„zn  (4.8) 

where  an  0  lie  in  the  circle 

\z\  <  1  4-  max  |at/fl„|  for  k  =  0, 1,  •  ,  n-1  (4.9) 

Applying  the  lemma  to  (4.3),  we  find  that  if  (u/,z)  is  a  zero  of  p(w,z),  the  following 
inequality  must  hold, 

|u>i  <  1  +  max\plc(z)/psl{z)\  for  k  =  0, 1,  •  ,  M  -  1  (4.10) 

If  w(z)  becomes  unbounded  then  we  need  that  Pu{z)  tend  to  zero.  Since  p.\/(z)  has  at 
most  N  zeros,  the  theorem  is  proven.  □ 

The  second  theorem  concerns  itself  with  branch  points,  the  second  type  of  singu¬ 
larity  considered. 

Theorem  4.2  The  complex  number  z0  is  a  branch  point,  i.e.,  two  roots  of  p(w,z), 
w,{z)  and  t n;(z),  are  equal  at  z0  where  w^(z0)  =  ui}(z0 )  =  m0;  if  and  only  if  the  partial 
derivative  of  p(w,z)  with  respect  to  vu  is  zero  when  evaluated  at  (u/0,2o). 

Proof:  We  will  denote  the  partial  derivative  of  p(w,  z )  with  respect  to  w  as  p,riu\  z). 
A  well  known  result  from  the  study  of  polynomial  functions  of  a  single  variable  states 
that  a  polynomial  p{z)  has  a  multiple  zero  at  z0,  if  and  only  if  z0  is  also  a  zero  of  the 
derivative  of  p(z)  [32,.  This  result  can  be  applied  trivially  to  the  case  of  a  bivariate 
polynomial.  If  (4.3)  when  evaluated  at  zn  has  a  multiple  root  at  w() ,  then  from  the 
discussion  above,  the  polynomial 

Pi(z.i)  2p2(zo)u>n  +  ■  •  ■  +  Mpn(z„)w^f~l  (4.11) 
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must  also  be  zero.  However,  this  polynomial  is  just  pv(w0,z0).  □ 


Example:  Consider  again  the  polynomial  of  the  previous  example.  In  this  case  we  can 
express  the  branches  of  w(z)  in  explicit  form. 


According  to  Theorem  4.1.  at  each  zero  of  p?(z)  =  1  -  r.  at  least  one  branch  of  tt>(z)  must  go 
to  infinity.  Here  the  zero  occurs  at  z  =  1.  and  we  see  from  (4.13)  and  Figure  4.1  that  w^(zj 
becomes  infinite  as  expected. 

The  polynomial  p,rtw.z)  for  our  example  is  given  by 

Pu-lw.z)  =  (-1  +  z  -  z2)  +  2(1  -  z)u>  (4.14) 

Calculating  p,r(w.z)  at  the  branch  point  z  —  .7913.  w  =  2.01G47.  we  find  that  £„(«;.  2)  is 
approximately  zero,  as  predicted  by  Theorem  4.2. 

From  (4.12).  we  see  that  all  branch  points  must  satisfy  the  polynomial  equation 

0  =  -3  +  62  —  5z2  +  2z3  +  z*  (4.15) 

Thus  for  this  example,  at  least,  p(w.z)  can  only  have  at  most  4  branch  points.  We  will  later 
characterize  more  generally,  the  number  of  branch  points  a  given  polynomial  may  have. 

Near  an  ordinary  point,  the  functions  w,(z)  can  be  associated  with  analytic  func¬ 
tions.  i.e.,  functions  possessing  a  Taylor  series  expansion. 


Theorem  4.3  Near  an  ordinary  point  z  =  a,  the  M  values  of  the  algebraic  function 
vl'(z)  are  defined  by  M  convergent  series 

w  =  u/l()  +  m,i(z  -  a)  -f  u>,;(2  -  a)2  +  ■  ■  •  (i  =  1,  •  •• ,  M)  (4.16) 

where  the  numbers  wl{,  are  the  M  distinct  roots  of  p(w,a)  =  0. 

Proof:  This  theorem  is  a  straightforward  application  of  the  Implicit  Function 

Theorem  (33,  p.  109( , 


Lemma  4.2  Let  F[w.z)  be  a  function  of  two  complex  variables  which  is  analytic  in  a 
neighborhood  z  -  z„.  <  r.  w  —  w(,\  <  p  of  the  point  ( ud<, ,  20 ) ,  and  suppose  that 


F[W„ ,  2„)  -  0. 


dF{ IU„,  2r,J 


(4.17) 
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Then  there  are  neighborhoods  M{z0),  M{w0)  suck  that  the  equation  F(w,z )  =  0  has 
a  unique  root  w  =  «/(«)  in  M(w0)  for  any  given  z  €  M(z0).  Moreover,  the  function 
w  =  w(z)  is  single-valued  and  analytic  on  M(z0)  and  satisfies  the  condition  w{z0)  =  w0. 

To  use  this  theorem  we  first  note  that  p(w,  z)  is  an  analytic  function  of  w  and  z  for 
all  finite  w,  z.  Furthermore,  from  Theorem  4.2,  if  (u>,0,  z0 )  is  an  ordinary  point,  then 
P«-{w ot  2o)  5^  0;  thus,  the  condition  of  (4.17)  is  satisfied.  We  conclude  then  that  a  locally 
analytic  unique  function  w(z)  exists  such  that  p(wi{z),z)  =  0  in  a  neighborhood  of  20 
and  such  that  u;,(20)  =  u>,o-  Since  there  are  M  solutions  to  the  equation  p{w,zQ)  =  0, 
there  are  M  such  functions.  □ 

A  corollary  of  this  theorem  is  that  each  u;,(z)  is  differentiable.  It  is  straightforward 
to  calculate  the  derivative  of  each  w,(z).  From  the  relationship, 


we  get  that 


Using  the  chain  rule, 


p(wi(z),z)  =  0 


dp(wi{z),z) 

dz 


p:(w,(z),  z)  +  pv(w,(z),  z)^^l  =  0 


(4.18) 


(4.19) 


(4.20) 


or 

dw,{z)  _  p;{wt{z),z) 

dz  pv{Wi{z),z)  (  ‘  ‘ 

The  crucial  importance  of  the  study  of  w(z)  is  that  each  branch  of  t u(z)  can  be 
associated  with  an  irreducible  factor  of  p(w,z).  This  relationship  is  stated  in  the 
following  theorem, 
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Theorem  4.4  For  each  branch  u>,-(z )  of  w(z)  given  by  (4-16)  and  each  ordinary  point 
z0  there  is  an  irreducible  factor  ofp(w,z),  p,(tv,z )  such  that 


I”  p,(u>,(z),  z)  =  0  for  all  z  in  a  neighborhood  of  z0  (4.22) 

Moreover ,  this  is  the  only  irreducible  polynomial  which  satisfies  the  equation  above. 

Proof:  Recall  that  any  polynomial  can  be  expressed  as  an  essentially  unique  prod¬ 
uct  of  irreducible  polynomials, 

p(w,  z)  =  px(w,  z)p2{w,  z)  ■  ■  •  pk(w,  z )  (4.23) 

Let  w,a  =  u/,(z0).  Since  p(wt0,Zo)  =  0,  we  have  according  to  (4.23), 

p(t";o.  z0)  =  Pi(«>,o,  z0)P2(u>.o,  z0)  •  •  •  Pt(u>.o,  z0)  =  0  (4.24) 

Thus  at  least  one  of  the  terms  p,(tofo,20)  must  be  zero.  However,  since  z0  is  an  ordinary 
point,  only  one  of  these  terms  can  be  zero.  Without  loss  of  generality,  let  pi(tt>,0,  z0)  =  0. 
Since  z0  is  an  ordinary  point  of  the  algebraic  function  of  p(w,z),  it  must  also  be  an 
ordinary  point  of  the  algebraic  function  corresponding  to  pl(w,z).  Thus  from  Theorem 
4.3  we  know  that  there  is  a  locally  analytic  function  g(z)  such  that  p,(g(z),z)  =  0 
everywhere  in  a  neighborhood  of  z0  and  ^(z0)  =  u>io-  However,  if  py{g{z),z)  =  0 
then  p(g(z),z)  =  0.  Since  according  to  Lemma  4.2,  the  function  tu,(z)  is  unique  in  a 
neighborhood  of  ztu  g{z)  and  w,(z)  must  be  the  same  in  this  neighborhood.  Therefore, 
p^w^z),  z)  =  0  for  all  z  in  a  neighborhood  of  z0.  That  this  is  the  only  such  irreducible 
polynomial  follows  from  the  fact  that  two  irreducible  polynomials  can  have  only  a  finite 
number  of  common  zeros.  A  theorem  to  this  effect  is  introduced  later  in  the  discussion 


of  Bezout’s  theorem. 


From  the  theorem  above,  we  see  that  we  can  associate  with  each  branch  of  t v(z) 
an  irreducible  factor  of  p(w,  z).  The  analogous  procedure  for  univariate  polynomials  is 
the  fact  that  we  can  associate  with  each  zero  a,  of  p(z),  an  irreducible  factor  of  p(z), 
namely  (a  -  a,).  The  crucial  difference  is  that  we  are  not  able  to  break  up  the  zeros 
of  p(tu,z)  into  isolated  zeros;  the  best  one  can  do  is  to  break  up  the  zeros  of  p(tn,a) 
into  different  continuous  sets  of  zeros.  The  reader  should  be  aware  that  each  branch  of 
w(z )  is  not  associated  with  a  distinct  factor  of  p(w,z),  that  is,  the  fact  that  w('z)  has 
M  branches  does  not  imply  that  p(w,z)  has  M  irreducible  factors.  For  example,  we 
could  have  that  both  Wi(z)  and  w^iz)  correspond  to  Pi(z). 

Although  we  will  have  the  opportunity  to  discuss  algebraic  functions  further,  the 
properties  described  above  are  sufficient  to  understand  the  rest  of  this  chapter.  In 
Chapter  5,  we  will  take  up  the  topic  again,  in  order  to  develop  some  more  results 
necessary  for  the  establishment  of  our  new  phase  retrieval  algorithm. 


4.3  Bivariate  Polynomial  Factorization 

4.3.1  Kronecker’s  Algorithm 

Although  univariate  polynomial  factorization  has  a  long  and  productive  history,  the 

factorization  of  bivariate  polynomials  has  enjoyed  relatively  little  progress.  The  first 

algorithm  for  factoring  a  polynomial  in  several  variables  is  due  to  Kronecker  in  1882 

34'.  Of  course,  Kronecker’s  interest  was  in  the  mathematical  problem  of  factoring 

polynomials  in  several  problems,  not  phase  retrieval.  The  basic  idea  is  to  convert  the 

bivariate  polynomial  into  a  univariate  polynomial,  and  then  factors  of  the  resulting 
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univariate  polynomial  are  associated  with  factors  of  the  original  bivariate  polynomial. 

Consider  a  general  polynomial  p(w,z)  of  degree  (M,N).  We  seek  to  find  a  poly¬ 
nomial  g(w,z )  which  is  a  factor  of  p(w,  z).  The  algorithm  consists  of  four  steps  as 
delineated  below: 

1.  Compute  degree  bound  Obtain  an  integer  d  >  max(M,  N). 

2.  Reduction  of  bivariate  to  univariate  polynomial  Generate  p{w)  which  is  given 

by 

p(w )  —  p(w,wd)  (4-25) 

Note  that  because  of  the  way  we  have  picked  d,  this  transformation  between  p(w) 
and  p{u),z)  is  reversible. 

3.  Factorization  of  univariate  polynomial  Factor  p(w)  into  irreducible,  i.e.,  lin¬ 

ear  factors, 

PM  =  ffiMfoM  ••  •  £,(tu)  (4.26) 

4.  Inverse  reduction  and  trial  division  From  the  factorization  in  step  3,  generate 

all  polynomials  which  divide  p(w).  Call  a  typical  factor  g(w).  Convert  g(vu)  into 
a  bivariate  polynomial  g{z,y )  via  the  relation 

g(w,wd)  =  g(w)  (4.27) 

and  check  if  g(w.z)  divides  p(w,z).  If  so,  stop;  if  not,  try  a  different  factor  If 

there  is  no  appropriate  gi w],  then  p(w,z)  is  irreducible. 

The  relationship  between  the  Canterakis  algorithm  and  Kronecker’s  algorithm  be¬ 
comes  apparent  if  we  note  that  the  polynomials  associated  with  fjnj  and  •■jm,  nj  in  the 


Canterakis  algorithm  are  related  by 

Pi  (tv)  =  p,(w,  w2X+l)  (4.28) 

Thus,  the  Canterakis  algorithm  picks  d  =  2 N  +  1  as  the  parameter  in  Kronecker’s 
algorithm.  Since  the  degree  of  pr(w,z)  is  (2N,2N)  for  an  N  by  N  image,  this  is  an 
appropriate  value  for  d.  Then,  the  algorithm  proceeds  to  factor  p.(ti>,  z )  via  Kronecker’s 
algorithm. 

A  severe  problem  with  Kronecker’s  algorithm  which  precludes  its  general  use  is  the 
tremendous  work  required  in  step  4  of  the  algorithm.  The  amount  of  computation 
increases  exponentially  as  the  size  of  the  polynomial  increases.  Since  the  Canterakis  al¬ 
gorithm  is  basically  Kronecker’s  algorithm,  it  also  suffers  from  this  exponential  growth. 

4.3.2  Kaltofen’s  Algorithm 

In  1982,  Kaltofen  [35],  presented  an  algorithm  for  reducing  a  bivariate  polynomial 
factorization  problem  to  a  univariate  factorization  problem  which  avoids  the  exponen¬ 
tial  growth  experienced  by  Kronecker’s  algorithm.  Again,  this  algorithm  was  developed 
solely  in  the  context  of  the  mathematical  problem  of  polynomial  factorization. 

The  first  step  in  Kaltofen’s  algorithm  is  to  find  an  ordinary  point  of  the  algebraic 
function  w(z).  It  can  be  supposed  without  loss  of  generality  that  20  =  0;  otherwise,  we 
just  need  to  make  a  linear  change  of  variables  in  z.  Thus,  it  can  be  assumed  that  a  xvu 
such  that  p(u>0, 0)  =  0  has  been  calculated  and  that  (tn0,0)  is  an  ordinary  point. 

The  second  step  of  the  algorithm  is  based  on  developing  a  power  series  expansion 
of  a  branch  of  the  algebraic  function  w(z)  around  zero.  That  such  an  expansion  exists 


is  guaranteed  by  Theorem  4.3, 


p(u>0  +  Viz  +  w2z2  +  - ,  z)  —  0  (4.29) 

There  are  many  ways  of  calculating  the  {in,}  set.  The  simplest  way  conceptually  is 
by  matching  powers  of  z,  i.e.  expanding  (4.29)  to  a  power  series  in  z  and  then  setting 
each  coefficient  in  the  power  series  to  zero.  This  yields  a  manner  of  evaluating  {w,} 
successively.  There  are  much  faster  algorithms  for  finding  the  {tn,}  set  which  are  based 
on  an  algebraic  version  of  Newton’s  method  [36,37];  however,  it  is  really  not  pertinent 
to  our  discussion  that  we  describe  these  methods  here. 

The  fourth  and  last  step  of  Kaltofen’s  algorithm  is  based  on  the  fact  that  each 
branch  of  t v(z)  is  characteristic  of  an  irreducible  factor  of  p(tn,  z)  by  virtue  of  Theorem 
4.4,  i.e.,  there  is  an  irreducible  polynomial  g(w,z)  that  divides  p{w,z)  and,  again  for 
all  a  in  a  neighborhood  of  z  =  0, 

g{w„  +  WiZ  +  w2z 2  +  ••  • ,  z)  =  0  (4.30) 

From  the  discussion  concerning  Theorem  4.4,  we  see  that  which  of  the  irreducible 
factors  of  p( w,  z)  is  specified  is  determined  by  which  irreducible  factor  of  p(w ,  z)  satisfies 
g(w„,  0)  =  0.  Picking  a  different  solution  tn„  may  result  in  a  different  irreducible  factor 
being  identified. 

Now  that  Taylor  series  defined  by  the  {u, }  set  has  been  extracted  and  associated 
this  series  with  one  of  the  irreducible  factors  of  g(w.z),  the  question  becomes  how  to 
extract  the  coefficients  of  g(tc,z )  given  the  coefficient  set  {in,}.  The  answer  is  to  use 
(4.30)  in  “reverse”,  i.e.,  calculate  the  coefficients  of  g{w,z )  from  the  known  {in,}  by 
matching  coefficients  in  z  and  setting  them  to  zero.  This  yields  a  set  of  homogeneous 


linear  equations  in  the  unknown  coefficients  of  g(w,z).  The  algorithm  tries  all  possible 
total  degrees  for  g{tu,  z)  up  to  totdeg(p)  until  it  finds  that  the  homogeneous  linear 
equations  have  a  non-trivial  solution.  In  the  reconstruction  from  Fourier  transform 
magnitude  problem  of  a  signal  with  an  irreducible  associated  polynomial,  we  would 
not  have  to  try  all  possible  sizes,  since  the  support  of  the  signal,  and  therefore  the 
degree  of  the  polynomial  factor,  is  known. 

Note  that  the  Kaltofen  algorithm  does  not  have  to  go  through  the  “combinatorial 
explosion”  required  by  the  Kronecker  algorithm.  Thus,  it  is  a  promising  vehicle  for 
factoring  polynomials  of  reasonable  size.  From  this  discussion  the  Kaltofen  algorithm 
seems  like  a  good  candidate  to  effect  the  factorization  of  the  polynomial  associated  with 
the  autocorrelation  function  as  was  done  in  the  Canterakis  algorithm,  but  without  the 
enormous  computational  requirements. 

Indeed,  we  have  tried  using  this  algorithm  and  have  found  that  it  successfully  re¬ 
constructs  signals  with  support  up  to  [0,  5]\  However,  as  is  well  known,  Taylor  series 
coefficients  have  to  be  calculated  exceedingly  accurately;  otherwise,  small  errors  in  the 
coefficients  lead  to  great  changes  in  the  function  as  the  series  is  evaluated  away  from 
the  origin.  Similarly,  the  Taylor  series  coefficients  in  Kaltofen’s  algorithm  have  to  be 
calculated  very  precisely,  or  else  the  equations  to  be  solved  do  not  match  the  actual 
solution.  We  have  found  that  these  round-off  effects  preclude  the  use  of  the  algorithm 
for  signals  with  support  larger  than  [0, 5j2.  Another  way  of  viewing  the  source  of  insta¬ 
bility  in  the  Kaltofen  algorithm  is  that  it  essentially  uses  local  information  to  extract 
the  factors  of  a  polynomial;  namely,  the  value  and  derivatives  of  the  locus  of  roots  of 
the  irreducible  polynomial  to  be  extracted. 


Chapter  5 


New  Closed  Form  Algorithm  for 
Reconstruction  from  Fourier 
Transform  Magnitude 


We  have  hinted  in  the  introduction  that  the  approach  we  consider  in  this  the¬ 
sis  is  to  explicitly  view  phase  retrieval  as  a  polynomial  factorization  problem  and  to 
search  for  an  algorithm  to  factor  polynomials  efficiently  and  accurately.  We  have  al¬ 
ready  noted  that  the  Canterakis  approach  for  reconstruction  from  Fourier  transform 
magnitude  is  implicitly  the  use  of  Kronecker’s  algorithm  to  factor  the  polynomial  asso¬ 
ciated  with  the  signal  autocorrelation  function.  However,  we  found  that  this  algorithm 
is  extremely  time  consuming  for  even  small  signals.  An  alternative  is  Kaltofen’s  al¬ 
gorithm  for  factoring  polynomials.  This  algorithm,  although  not  as  computationally 
expensive  as  Kronecker’s  algorithm,  suffers  from  severe  numerical  instability  due  to  the 
use  of  information  which  is  very  sensitive  to  round-off  effects. 

In  this  chapter  a  new  algorithm  for  factoring  bivariate  polynomials  is  described. 
Although  the  algorithm  is  applicable  to  the  factorization  of  any  bivariate  polynomial, 
it  is  specifically  developed  as  a  vehicle  for  the  phase  retrieval  problem.  This  algorithm 
does  not  increase  in  complexity  as  the  region  of  support  increases  as  quickly  as  the 


Canterakis  or  Deighton  algorithms  do  and  is  more  numerically  stable  than  applying 
the  Kaltofen  algorithm,  allowing  us  to  factor  larger  polynomials  and  thus  reconstruct 
larger  images  from  the  Fourier  transform  magnitude. 

The  algorithm  is  based  on  analyzing  the  set  of  complex  pairs  (w,  z)  where  the 
polynomial  associated  with  the  autocorrelation  function  is  zero.  Although  the  use  of 
this  set  of  zeros  to  factor  pr[w,z )  is  original,  the  importance  of  this  set  of  zeros  to  the 
problem  of  phase  retrieval  has  long  been  recognized,  beginning  with  the  early  work  of 
Napier  and  Bates  [38] . 

5.1  Background  and  Overview 

Before  developing  the  algorithm  formally,  it  is  helpful  to  sketch  the  general 
idea  behind  the  method.  Consider  a  reducible  polynomial  p{w,z )  with  two  irreducible 
factors  r(w,z )  and  s(w,z );  thus 

p(w,  z)  =  r(w,  z)s(w,  z)  (5.1) 

If  (u),  z)  is  a  zero  of  s{xv,  z)  then  it  must  also  be  a  zero  of  p(tn,  z).  However,  if  ( w ,  z )  is 
a  zero  of  p(w,z)  then  a  priori  one  cannot  tell  whether  this  zero  corresponds  to  r(w,z) 
or  s(w,  z ).  Of  course,  (w,  z)  must  be  a  zero  of  at  least  one  of  them.  Generating  zeros  of 
p(tu,z)  is  easy  since  p(w,z)  is  known;  the  difficulty  resides  in  assigning  each  calculated 
zero  of  p(w,z)  to  the  appropriate  irreducible  factor  of  p(w,z ).  The  main  objective  of 
the  development  discussed  below  is  to  generate  such  an  assignment  strategy,  so  that 
eventually  we  have  a  large  number  of  zeros  which  are  all  guaranteed  to  correspond 
exclusively  to  either  r(w,z)  or  s(w,z). 

Suppose,  then,  that  we  nave  at  our  disposal  a  large  number  of  zeros  of,  for  example, 
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s(tb',z).  Then,  for  each  zero  (wk)zk),  the  following  equation  must  hold, 

s{wk,zk)  =  0  =  Soo  +  Sl0wk  s01zk  +  snwkzk  H -  (5.2) 

where  si;  are  the  coefficients  of  s(w,  z ),  This  equation  is  linear  in  the  unknown  s,y.  By 
examining  the  null  vectors  of  this  set  of  equations  we  will  be  able  to  reconstruct  the 
■coefficients  of  s(w,z )  and  thus  s(w,z )  itself. 

5.1.1  New  Factorization  Algorithm 

An  interpretation  of  w(z),  the  algebraic  function  corresponding  to  p(w,z),  which 
will  prove  to  be  crucial  in  our  later  development,  is  to  consider  it  as  a  multi-valued 
mapping  of  a  path  in  the  z-plane  to  several  paths  in  the  tn-plane.  A  path  is  a  complex¬ 
valued  function  of  a  real  parameter  t,  where  t  varies  from  t0  to  tv 

z(t)  =  x(t)+jy(t)  (5.3) 

such  that  x{t)  and  y{t)  are  continuous.  Since  we  will  be  free  to  choose  2(1),  it  will  be 
convenient  to  assume  also  that  x(t )  and  y(t)  are  also  piece-wise  differentiable  functions 
of  t.  The  mechanism  by  which  this  mapping  occurs  is  easy  to  visualize.  Consider  an 
ordinary  point  z„  =  z( 0)  and  all  M  solutions  to  the  (univariate)  polynomial  equation 
p(tL\zt ,)  =  0.  We  denote  each  of  these  solutions  by  u/,(0)  for  t  =  1,  ■ ,  M,  Figure 

5.1.  As  z„  is  changed  in  a  continuous  manner,  thus  generating  z(t),  the  corresponding 
roots  of  the  equation  p(w,z{t))  will  also  change  in  a  continuous  manner,  resulting  in 
the  paths  Figure  5.2  This  view  of  w{z 1  mapping  paths  in  the  2-plane  to  paths 

in  the  u-nlane  is  supported  by  the  following  theorem  10,  p.25] , 


Theorem  5.1  If  z  =  z(t)  (ft  <  i  <  t2)  is  a  path,  consisting  entirely  of  ordinary  points 
of  an  algebraic  function  w(z),  then  the  values  of  t v(z)  along  z(t)  form  a  set  wk(t ) 
(<!  <  t  <  t2;  k  —  1,  • ,  M)  of  M  paths. 

In  terms  of  the  M  branches  of  w(z),  the  paths  described  in  this  theorem  are  given 

by, 

t i>k{t)  =  Wk{z{t ))  (5.4) 

Thus,  the  trajectory  of  each  of  the  paths  is  specified  by  each  of  the  M  branches  of 
ru(z). 

Example:  As  an  illustration,  consider  a  path  in  the  2-plane  given  by  z(t)  =  ej2r,.(0  <  t  < 
1).  Figure  5.3.  The  polynomial  under  study  is  p(w.z)  =  to2  -  z.  We  now  evaluate  the  two 
paths  tei(t)  and  u)o(t)  which  form  the  image  of  z(t)  under  w(z).  At  each  value  of  t.  we  have, 
p(w,(t).  z(t))  =  0,  or  wi(t)  =  e>Kt.  wi(t)  =  e~’*‘ .  The  resulting  paths  are  displayed  in  Figure 
5.4. 

Theorem  5.1  extends  our  notion  of  t v(z)  from  the  local  characterization  given  by 
Theorem  4.3  to  a  global  one.  Note  that  this  global  characterization  comes  at  the 
expense  of  no  longer  being  able  to  consider  each  u>,(a)  individually. 

Example :  Consider  the  two  paths  za(t )  and  zi,(t)  (0  <  t  <  1)  given  by 

z0(t)  =  e’2*‘  (5.5) 

zb(t)  =  2  -  e.’ut  (5.G) 

They  are  depicted  in  Figure  5,5.  We  consider  again,  the  polynomial  p(w.z)  =  u-2  —  z.  Since 

2„(0)  =  2j(0).  we  can  associate  with  za{t )  and  zb(t)  two  paths  which  begin  at  the  same  value 

of  w  =  l.  w„(t)  and  u’(,(t).  The  corresponding  paths  are  drawn  in  Figure  5.G.  Although 

2„(1)  =  zb(  1 ) ,  we  find  that  ten(l)  ^  wj(l).  i.e..  we  have  “moved"  from  one  branch  of  w(z)  to 

another. 

The  reader  may  at  this  point  realize  the  similarity  between  the  Figures  5.3  to  5.6 
and  the  technique  of  root  locus  analysis  which  forms  an  important  part  of  classical 
control  analysis.  The  observation  is  well-justified  since  root  locus  analysis  seeks  to  find 
the  zeros  of  the  function 

P{s,K)  =  H{s)  +  KG{s)  (5.7) 
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where  H(s)  is  the  feed-forward  transfer  function  numerator,  G(s)  is  the  feed-forward 
transfer  function  denominator  and  K  is  the  feedback  gain.  The  similarity  becomes 
obvious  once  we  consider  (5.7)  as  a  polynomial  in  the  two  variables  s  and  K.  Moreover, 
it  is  easy  to  show  that  P(s,  K)  is  an  irreducible  polynomial  as  long  as  H(s)  and  G(s) 
do  not  have  a  common  factor.  Thus,  many  of  the  observations  we  have  made  so  far  and 
will  be  making  in  the  future  also  pertain  to  root-locus  analysis.  Where  our  discussion 
and  root-locus  analysis  depart  is  that  although  in  root-locus  analysis,  K  is  restricted 
to  be  real,  we  will  allow  both  w  and  z  to  take  on  complex  values. 

We  need  one  more  result  regarding  the  path  mapping  property  of  algebraic  func¬ 
tions.  This  next  result  is  analogous  to  Theorem  4.4  which  was  used  by  Kaltofen  to 
develop  his  bivariate  polynomial  factorization  algorithm.  Recall  that  Theorem  4.4 
stated  simply  that  each  branch  t u,[z)  of  the  algebraic  function  t u(z)  corresponding  to 
a  polynomial  p(u>,  z)  can  be  locally  associated  with  an  irreducible  factor  of  p{w,  a), 
i.e.,  there  is  an  irreducible  polynomial  p,(u;,  z)  which  divides  p(xv,z)  and  for  which 
p,(w,(z),  z)  =  0  is  true  everywhere  in  a  neighborhood  of  an  ordinary  point.  This  next 
theorem  states  a  similar  result  for  each  path  u/,-(t)  which  is  an  image  of  z(t)  under 
w(z). 


Theorem  5.2  Let  z{t)  be  a  path  in  the  z -plane  consisting  exclusively  of  ordinary  points 
of  an  algebraic  function  w(z)  corresponding  to  a  polynomial  p(w,z).  If  the  path  i e\(f)  is 
an  image  of  z(t)  via  w(z )  in  a  given  interval,  then  there  is  an  irreducible  factor  p,(w,z) 
which  divides  p(u>,z)  such  that  p,(u;,(£),  z(t))  =  0  for  all  t  in  the  given  interval. 


Proof:  This  follows  rather  directly  from  the  fact  that  z(t)  only  includes  ordinary 
points  and  that  z(t)  and  all  the  polynomials  involved  are  continuous  functions.  Suppose 
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that  for  t  =  t0,  (w, (to),z{to))  is  a  zero  of  pi(w,z)  where  p,(tn, z }  is  an  irreducible  factor 
of  p(u>,  z).  Express  p(w ,  z)  as  the  product  p,( w,  z)q( w,  z).  We  want  to  show  that  for  all 
subsequent  values  of  t,  (u>,(£),  z(£))  is  a  zero  of  p,(u/,  z).  We  will  have  to  make  use  of 
Theorem  5.1,  i.e.,  if  z(£)  consists  of  only  ordinary  points,  then  tn(£)  is  also  a  continuous 
path. 

Suppose  that  at  some  point,  the  hypothesis  is  not  true,  i.e.,  p,(w,(t),  z(t))  /  0,  say 
at  <i.  Since  (in,(£),  z(t))  is  a  zero  of  p(to,  z)  for  all,  we  must  have  that  g(m(t1),z(t1))  = 
0.  Because  of  continuity,  there  must  be  a  transition  t'  between  t0  and  £x  such  that 
p,(w(t’),  z(t'))  =  q(w(t'),  z(t'))  =0.  However,  this  would  mean  that  pv(w(t" ),  z(t*))  = 
0  or  that  z(V)  is  a  singular  point  of  w(z).  This  contradicts  the  assumption  that  z(t) 
consists  of  only  ordinary  points.  The  fact  that  this  irreducible  polynomial  p,(w,z)  is 
the  only  one  which  contains  (i u,(t),  z(t))  as  a  zero  for  all  t  will  become  obvious  once  we 
introduce  Bezout’s  theorem  later  on.  □ 

The  theorem  above  is  exactly  what  is  needed.  From  Theorem  5.2  we  see  that  the 
zeros  of  p(u»,z)  which  lie  on  a  path  pair  (u>(t),z(£))  must  all  be  from  one  irreducible 
factor  of  p(w,z).  An  algorithm  for  isolating  an  irreducible  factor  of  p{w,z)  can  be 
described  in  the  following  three  steps: 

•  Pick  a  path  z(t )  consisting  entirely  of  ordinary  points  of  p{w,z). 

•  Find  a  path  w(t)  such  that  p(w(t),z(t))  =  0  for  all  t  of  interest. 

•  Find  the  irreducible  polynomial  d(w,z)  such  that  d(w(t),  z(t))  =  0  for  all  t  of 
interest.  According  to  Theorem  5.2,  the  resulting  polynomial  d{w,z)  will  divide 
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5.1.2  Picking  a  Path  Consisting  of  Ordinary  Points 

In  order  to  study  which  paths  in  the  z-plane  consist  of  only  ordinary  points,  we  need 
to  discuss  in  more  detail  the  behavior  and  propensity  of  points  in  the  z-plane  which  are 
not  ordinary,  i.e.,  the  singular  points.  Recall  that  there  are  two  types  of  singular  points: 
the  first  type  is  composed  of  all  z  such  that  one  branch  tn,(z)  becomes  infinite.  We 
have  already  shown  in  Theorem  4.1  that  there  are  degv(p )  such  z.  The  second  type  of 
singular  point  consists  of  all  z  such  that  two  branches  coalesce  at  z,  i.e.,  w,(z)  =  Wj(z) 
at  z.  Theorem  4.2  stated  that  if  (tt>o,2o)  *s  such  that  w0  =  ti>,(z0)  =  Wj(z0),  then  not 
only  p(w0,z0 )  =  0  but  also  pv{w0,  z0)  =  0  as  well.  Since  p„(tu,  z)  is  also  a  bivariate 
polynomial  and  since  (u>0,Zo)  resides  in  the  common  zeros  of  the  two  polynomials 
p(xv,  z)  and  p„.(t/;,  z),  we  need  to  discuss  briefly  the  topic  of  the  intersection  of  the  zeros 
of  bivariate  polynomials. 

In  general,  two  bivariate  polynomials  will  only  have  a  finite  number  of  common 
zeros.  This  can  be  supported  intuitively  by  the  fact  that  a  possible  common  zero 
(u;0,z0)  of  p(u>,z)  and  q{w,z)  can  be  represented  by  the  four  dimensional  real  vector 
[Re(u’0),  7m(u;()),  J?e(z0),  7m(z0)]T.  On  the  other  hand,  (u>0,Zo)  must  satisfy  simultane¬ 
ously  the  four  equations  below, 

Re(p(w0,  20))  =  Im(p(w0,z0))  =  Re(q(w0,  z0))  =  Im(q(w0,z0))  =  0  (5.8) 

The  maximum  number  of  such  common  zeros  is  the  subject  of  Bezout’s  theorem  below 
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Theorem  5.3  If  p{w,z )  and  q(w,z)  are  bivariate  polynomials  of  total  degree  r  and 
s  respectively  with  no  common  factors,  then  there  are  at  most  rs  distinct  pairs  (tn,  z) 
where 


p{w,z)  =  q{w,z )  =  0 


(5.9) 


Example:  All  example  where  the  maximum  number  of  common  zeros  is  attained  with  finite 
w  and  z  is  the  pair  of  polynomials 


p(w.  z)  =  w2  +  z2  -  4  (5.10) 

q(w.z)  =  w2  -  z2  -  1  (5.11) 

The  set  of  zeros  of  these  polynomials  for  real  (u>.z)  is  plotted  in  Figure  5.7.  The  total  degree 
for  each  polynomial  is  2.  The  common  zeros  occur  at  (±y/3.  ±\/5).  Since  the  total  degree  of 
each  polynomial  is  2.  Bezout's  theorem  predicts  a  maximum  number  of  common  zeros  of  4. 
which  is  attained  exactly  in  this  case. 

One  application  of  Bezout’s  theorem  which  we  will  have  need  to  make  use  of  later 
on,  is  that  it  specifies  the  minimum  number  of  zeros  needed  to  specify  an  irreducible 
polynomial  to  a  scale  factor.  We  will  discuss  this  point  in  more  detail  in  a  later  part 
of  this  chapter. 

Bezout’s  theorem  applies  to  any  two  polynomials  with  the  specified  total  degree. 
In  much  of  our  discussion  we  will  know  the  degree  in  each  variable  as  well  as  the 
total  degree  of  the  polynomial  in  question.  An  obvious  question  then  is  whether  we 
can  somehow  tighten  the  bound  on  the  number  of  common  zeros  if  the  degree  of  the 
polynomials  involved  is  known.  In  a  sense  we  cannot  do  this  because  there  is  a  strong 
form  of  Bezout’s  theorem  where  it  is  shown  that  the  bound  described  above  is  actually 
achieved.  If  the  degree  of  p(w,  z)  is  ( MP,NP )  with  the  z-Vp  term  non-zero  and  the 
degree  of  q(vj,z)  is  (A/,,  Nq)  with  the  w*l*zx,>  term  non-zero  then  the  strong  form  of 
Bezout’s  theorem  states  that  there  are  ( Mp  +  NP){MI}  4-  IV,)  common  zeros.  However, 
many  of  these  zeros  are  “phantom”  zeros  where  w  or  z  or  both  are  infinite.  In  our 
study  we  will  be  concerned  with  finite  zeros.  Therefore,  it  is  important  to  be  able  to 
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study  the  number  of  finite  zeros  which  two  polynomials  with  specified  degree  can  have 
in  common.  This  is  the  topic  of  Appendix  B  where  the  following  theorem  is  proven, 


Theorem  5.4  Let  p(w ,  z )  and  q{w,  z)  be  two  polynomials  of  degree  (Mp,  Np)  and  ( Mq ,  Nq) 
with  no  common  factors,  then  there  are  at  most  NpMq  4-  NqMp  pairs  of  finite  complex 
numbers  ( w,z )  such  that 

p(w,z)  =  q(w,z)  =  0  (5.12) 


That  the  bound  above  can  be  attained  exactly  for  some  p(w,z)  and  q(w,z)  is 
demonstrated  by  the  following  example: 

Example:  Consider  the  two  polynomials  of  degree  (1. 1) 

p(w,  z)  =  wz  -  2  (513) 

q(w,  z)  =  wz  +  w  —  z  -  3  (5.14) 

The  zero  sets  for  w.  z  real  are  drawn  in  Figure  5.8.  Since  these  two  polynomials  have  2  common 
finite  zeros,  the  bound  of  Theorem  5.4  is  achieved. 

Recall  that  branch  points  of  p(w,  z)  reside  in  the  intersection  of  the  zeros  of  p(w,  z) 
and  pv(w,z).  The  theorems  just  discussed  specify  the  number  of  zeros  in  common 
between  two  polynomials  which  are  relatively  prime.  As  soon  as  we  show  that  in  most 
cases  p(w,  z )  and  p,r(,  z )  are  relatively  prime,  then  we  can  use  these  theorems  to  specify 
the  maximum  number  of  branch  points  p(w,z)  may  have,  which  in  turn  bounds  the 
number  of  singular  points  in  the  z-plane. 


Theorem  5.5  Letp(w,  z)  be  expressed  as  a  product  of  irreducible  primitive  polynomials 
and  a  constant 

r 

p{w,z)  =  a  U  p,{w,z)  (5.15) 

i=i 

If  p,(w,z)  -i  pa.(w,z)  for  i  =£  k,  then  p(w,z)  and  pv{w,z )  are  relatively  prime. 
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Figure  5.8:  Intersection  of  zeros  of  t vz  -  2  and  wz  ~  w  —  z  -  3. 
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Proof:  A  polynomial  such  that  its  irreducible  factors  dc  not  repeat  is  called  a 
square-free  polynomial.  Thus,  the  condition  above  is  that  p(w,  z)  be  square-free.  As¬ 
sume  that  p(w,  z)  and  pu.(u»,2)  have  a  common  non-trivial  irreducible  factor  Pi(w,z). 
Therefore  pi(w,z)  divides  pv{ u>,z),  say  pte(w,z)  =  pl(tu,z)q{w,z).  However  pv(w,z) 
can  be  expressed  as, 


P i..r(u>>  z)A(w,  z)  =  (g(tu,  2)  -  B{w,  2)]  p^w,  z)  (5.18) 


Therefore  pt( w,  z)  divides  A{w,  z )  or  pi.a.(tu,  2).  However,  since  p(tt>,  2)  has  no  repeated 
factors,  A(w,z)  and  p^tn,  z)  are  relatively  prime.  Furthermore,  p,.,r(u>,2)  is  of  lower 
total  degree  than  pt  ( w ,  2)  and  therefore  cannot  contain  Pi(w,  z)  as  a  factor  since  pi(tn,  2) 
is  irreducible.  Therefore  pl.,r{xv,z)A{w,z)  and  p^w^z)  are  relatively  prime.  Thus  the 
assumption  that  p{xv,z)  and  pir(w,z)  have  a  common  factor  leads  to  a  contradiction. 


Using  Theorem  5.5,  we  can  deduce  that  for  bivariate  polynomials  which  are  square- 
free,  the  number  of  singular  points  in  the  2-plane  is  finite,  and  in  fact  the  maximum 
number  of  finite  singular  points  can  be  bounded  by  the  minimum  ol  the  bounds  of 
Theorems  5.3  and  5.4. 
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Theorem  5.6  Let  p(w,  z )  be  a  square-free  polynomial  of  degree  (M,  N)  and  total  degree 
npi  then  the  number  of  finite  singular  points  is  bounded  by 
min{np[np  —  1 ) ,  2 (M  -  l)N)  +  N. 

Proof:  Since  p(w,z)  is  square-free,  p(w,z)  and  pw(w,z )  are  relatively  prime.  Since 
the  total  degree  of  pu.(w,z)  is  np  —  1  and  the  degree  of  pv(w,  z)  is  ( M  —  1 ,  JV),  the 
maximum  number  of  finite  common  zeros  is  min(np(np  —  1 ) ,  2 (M  —  1  )N).  The  first 
term  is  due  to  Bezout’s  theorem;  the  second  term  is  due  to  our  intersection  Theorem 
5.4.  From  Theorem  4.2,  this  means  that  the  number  of  branch  points  of  p{w,z)  is 
bounded  by  min(np(np  —  1),  2 (M  -  l)iV).  From  Theorem  4.1,  the  number  of  singular 
points  where  a  branch  of  t v(z)  becomes  infinite  is  at  most  N.  Combining  the  two  sets 
of  singular  points  we  conclude  with  the  total  bound  of  min(np(np  —  1),  2MN  —  A^)  +  N. 

The  importance  of  the  theorem  above  at  this  point  is  not  the  actual  bound  on  the 
number  of  singular  points,  but  rather  that  the  number  is  finite.  This  means  that  if  we 
pick  a  path  z(t)  at  random,  it  will  almost  surely  consist  exclusively  of  ordinary  points. 

The  previous  discussion  showed  that  almost  all  paths  in  the  0-plane  consist  of  only 
ordinary  points,  since  singular  points  form  a  very  sparse  set,  consisting  of  only  a  finite 
number  of  points.  In  the  following  section  we  address  the  problem  of  specifying  a  path 
u;(t)  for  a  given  path  z(t)  such  that  p(w(t),  z(t))  is  zero  for  all  t  of  interest. 

5.1.3  Tracking  Zero  Paths 

Given  a  path  z(t)  consisting  of  ordinary  points,  and  an  initial  zero  of  p(w,z),  say 
p(c,  0(0)),  it  is  straightforward  to  pose  a  differential  equation  that  u>(t)  must  satisfy. 
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Recall  from  (4.21)  we  calculated  the  derivative  of  each  branch  of  w(z)  as, 

dwj(z)  _  p:(wj(z),z ) 
dz  p,r{wt[z),z) 

This  in  turn  implies  that  the  derivative  of  w(t)  with  respect  to  t  must  be 

dw(t)  _  p.(m(f),z(f))  dzjt) 

dt  pv(w(t),  z(t))  dt 


(5.19) 


(5.20) 


We  are  ieft  with  the  following  first  order  initial  value  problem: 
Find  the  function  w(t)  which  satisfies,  uj(0)  =  c  and 


dw(t)  _  p,(w(t),  z(t))  dz(t) 

dt  p.{tv(t),  z{t))  dt 


By  solving  the  initial  value  problem  above,  we  can  calculate  a  w(t)  for  any  z[t)  such 
that  p(w(t),  z(t))  =  0  for  all  t  of  interest. 


5.1.4  Extracting  Polynomial  Coefficients  from  Zeros 


At  this  point  we  have,  at  least  conceptually,  an  arbitrarily  large  number  of  ze¬ 
ros  of  one  irreducible  factor  of  p(w,z ),  say  d(w,z);  namely  all  the  complex  pairs 
( w(t),z(t )).  The  next  part  of  the  procedure  extracts  the  coefficients  of  d(w,  z)  from 
its  zeros.  This  problem  was  first  considered  and  solved  by  Curtis  [39,40] .  Basically, 
we  sample  {w(t),  z{t))  at  4,  t  =  1,  ••  ,  K  to  yield  K  simultaneous  homogeneous  linear 
equa'lons  for  the  coefficients  of  d(w,z), 


(irg-i  ('/)  'ifQttii | 


m=D  r»=0 


w™z%  =  0  for  k  =  1, 


where  d,n.„  is  the  set  of  coefficients  of  d[w,z). 
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There  are  two  points  that  need  to  De  addressed  at  this  point.  First  is  the  question 
of  how  many  samples  of  the  zeros  of  d(tv,z)  are  needed  to  guarantee  that  there  will  be 
essentially  only  one  solution  to  (5.22).  The  second  question  is  how  to  specify  degv(d) 
and  deg.(d)  without  knowing  d(w,z). 

The  question  of  determining  the  degree  of  d(w,  z)  is  addressed  in  Appendix  C.  The 
reason  why  we  do  not  discuss  the  problem  here  is  that  in  the  case  of  the  phase  retrieval 
problem,  the  issue  does  not  arise  as  we  will  see  later.  Let  us  assume  then  that  the 

degree  of  the  irreducible  factor  of  p{w,z)  to  be  extracted  is  known.  We  now  want  to 

specify  K  such  that  the  coefficients  of  d(w,z)  will  be  retrievable  from  (5.22). 

Consider  a  non-zero  solution  to  (5.22),  say  the  coefficients  am  „.  We  can  associate 
with  this  set  of  coefficients  a  polynomial  a(u>,z)  given  by, 

degu  (rf)  deg  Ad) 

a{w,z)  =  Yl  Yl  a.m.nwmzn  (5.23) 

m=0  n— 0 

Since  d(w,z )  is  irreducible,  and  a(w,  z)  and  d(w,z )  have  the  same  degree,  d(w,z) 
and  a(w,z )  must  be  relatively  prime  or  related  by  a  constant  factor.  From  Theorem 
5.4,  then  d(w,z)  and  a(w,z)  can  have  at  most  2degv(d)deg:(d)  zeros  in  common.  Since 
am  n  is  a  solution  to  (5.22),  the  polynomial  a(w,z )  must  satisfy, 

a{wk,zk)=0  (5.24) 

Thus,  afw.z 1  and  d{w,z )  have  K  zeros  in  common.  If  we  pick  K  >  2 degl{d)degv(d), 
then  a(w,z)  must  be  related  to  d{w,z)  by  a  constant  factor  which  implies  in  turn  that 
a,„  „  and  d„,  „  are  related  by  a  constant  factor.  Therefore,  by  finding  a  nuil  vector  of 
the  matrix  displayed  in  (5.22)  we  have  recovered  an  irreducibie  factor  of  p(w,z).  Note 


that  using  Theorem  5.3,  the  minimum  number  of  equations  required  to  guarantee  the 
correct  reconstruction  is  twice  the  number  specified  by  Theorem  5.4. 

Let  us  review  the  steps  necessary  to  find  an  irreducible  factor  of  a  given  square-free 
bivariate  polynomial  p(w,  z). 

1.  Find  an  ordinary  point  z0  of  the  algebraic  function  corresponding  to  p(to,z). 

2.  Using  z0 ,  find  a  w0  such  that  (w0,z0)  is  a  zero  of  p(w,  z).  Such  a  w0  can  be 
calculated  via  the  method  described  in  Section  4.2.2.  This  complex  pair  will  also 
be  a  zero  of  an  irreducible  factor  of  p(w,z),  d(w,z). 

3.  Generate  a  path  z(t)  which  begins  at  z0  and  which  consists  only  of  ordinary 
points.  Since  a  square-free  polynomial  has  a  finite  number  of  singular  points, 
almost  any  path  will  consist  of  only  ordinary  points. 

4.  Use  the  differential  equation  (5.21)  and  the  initial  condition  w0  =  u;(0)  to  calcu¬ 
late  a  path  w(t)  which  satisfies  p(w{t),  z(f))  =  0  for  all  t  of  interest. 

5.  Sample  the  path  pair  (to(£),  z(£))  to  yield  ( tu*,zfc )  and  solve  (5.22)  to  get  the 
coefficients  of  an  irreducible  factor  of  p(w,z),  d{xv,z). 

5.2  Considerations  for  Phase  Retrieval 

In  the  discussion  above  we  considered  the  polynomial  factorization  problem  in  gen¬ 
eral.  Here  we  would  like  to  consider  the  phase  retrieval  problem  specifically.  We  first 
note  that  since  in  almost  all  practical  cases,  the  polynomial  p,(u;,z)  is  irreducible,  the 
autocorrelation  polynomial  pr(w,z)  will  only  have  two  factors,  namely  pr(w,z)  and 


Pj(w,z).  Therefore,  unlike  the  general  polynomial  factorization  problem,  the  degree  of 
the  factors  is  known;  if  deg{pr )  =  (2M,2N)  then  deg{px)  =  deg(px)  =  (M,N).  One  fi¬ 
nal  point  to  consider  is  under  what  situations  p,(w,  z)  is  a  square-free  polynomial  which 
implies  in  turn  that  pr(w,z)  has  a  finite  number  of  singular  points.  Since  in  practical 
cases,  px(w,z)  is  irreducible,  we  restrict  our  attention  to  this  case  only.  In  this  situa¬ 
tion,  pr(w,z)  has  only  two  factors,  px(w,z)  and  px(u>,  z).  Thus,  if  px(w,z)  ^  px(w,z) 
then  pr{ii\z)  is  square-free.  Pictorially,  this  means  that  x[m,  n]  when  rotated  by  180 
degrees  does  not  look  identical  to  the  original  image.  Thus,  we  have  the  following 
theorem, 

Theorem  5.7  If  the  polynomial  associated  with  x[m,n\,  px(w,z )  is  irreducible,  and 
px(w,z)  is  not  equal  to  its  mirror  polynomial,  then  the  polynomial  associated  with 
r[m,  nl  will  be  square-free. 

Since  most  images  do  not  have  this  symmetry  property  we  can  be  assured  that 
p,(w,  z )  will  be  square-free  and  therefore  it  will  have  a  finite  number  of  singular  points. 
This  implies  in  turn  that  most  paths  in  the  z-plane  consists  of  only  ordinary  points. 

An  M  by  N  pixel  image  will  have  an  associated  polynomial  of  degree  (A/—  1 ,  JV-  1). 
Thus,  if  this  polynomial  is  irreducible,  it  is  specified,  according  to  Theorem  5.4,  by  a 
number  of  zeros  greater  than  2 (M  -  1  )(N  -  1).  The  minimum  number  of  equations  in 
(5.22)  necessary  to  assure  a  unique  solution  becomes  2 (M  —  1)(AT  —  1)  4-  1. 

The  total  algorithm  for  reconstruction  from  Fourier  transform  magnitude  of  a  signal 
with  an  irreducible  associated  polynomial  can  be  decomposed  into  the  following  steps. 


1.  Calculate  r’m,n]  from  ;X(e,u,  via  (2.13). 


2.  Form  pr( to,  2)  from  r[m,  n  . 


3.  Pick  z 0  such  that  it  is  an  ordinary  point  of  pf(tt;,z).  Almost  all  values  of  z()  will 
be  adequate. 

4.  Find  a  w0  such  that  pr{w0,z0)  =  0. 

5.  Pick  a  path  in  the  2-plane  consisting  of  ordinary  points  such  that  2(0)  =  z0  and 
solve  the  differential  equation  using  w0  as  the  initial  condition  and  (5.21). 

6.  Pick  enough  ( w ,  z )  pairs  such  that  (5.22)  has  a  unique  solution.  By  our  version 
of  Bezout’s  theorem  the  minimum  number  is  2 (M  —  1  )(N  —  1)  +  1  for  an  M  by 
N  image. 

7.  Since  p,(u>,  z)  is  irreducible,  then  the  solution  to  the  set  of  equations  above  yields 
coefficients  proportional  to  the  image  values  z[m,n)  or  x\-m,  -n]. 

8.  Normalize  the  coefficients  found  such  that 

^z[m,n]J  =  r[0,0)  (5.25) 

m  .n 

A  conceptual  flowchart  of  the  steps  involved  in  the  algorithm  is  depicted  in  Figure  5.9. 


Example:  Before  considering  some  practical  images,  it  is  instructive  to  illustrate  the 

algorithm  using  a  small  example.  Consider  the  situation  where  we  are  to  reconstruct  the 
following  imatre  from  the  Fourier  transform  magnitude. 

4  2  1 

x[m.  n]  =  0  3  3  (5.20 

3  0  1 


This  image  has  associated  polynomial  given  by. 

p,(tw. z)  =  4  +  iz  +  z2  +  3i vz  +  3u>22  +  3te2  4-  :2u'2  (5.27) 


FTM 


autocorrelation  function  =  r[m,  n] 


associated  polynomial  =  pr(w,z) 


zeros  of  pr(w,z ) 


o 


zeros  of  px(w,z)  zeros  of  px(w,z) 


m. 


—  m,  — 


As  a  path  in  the  2-plane.  we  will  use  z(t)  =  .&e}2xt ,  for  t  between  0  and  1.  i.e..  a  circle  of  radius 
.8  centered  at  the  origin.  Figure  5.10.  The  reason  for  using  this  path  is  discussed  in  the  next 
section.  For  this  path,  we  can  calculate  the  paths  u\(t)  in  the  w-plane  such  that  2(t))  = 

0  for  all  t  of  interest.  These  paths  are  drawn  in  Figure  5.11.  We  can  similarly  calculate  the 
paths  in  th  w-plane  which  correspond  to  zeros  of  the  mirror  polynomial  of  p,(ti>.z),  i.e.,  the 
polynomial  associated  with  the  signal. 

1  0  3 

3  3  0  (5.28) 

1  2  4 

These  paths  are  shown  in  Figure  5.12.  Comparing  Figures  5.il  and  5.12.  we  see  that  the  paths 
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Figure  5.10:  Path  chosen  in  the  2-plane  for  reconstruction  example.  The  unit  circle  is 
drawn  as  a  dotted  line. 
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Figure  5  11:  Resulting  paths  in  the  iy-plane  for  given  z(t).  The  path  u>i(t)  is  drawn  as 
a  soiid  line,  while  is  drawn  as  a  dashed  line. 


are  not  reciprocal  images  as  might  be  expected.  This  would  only  be  the  case  if  the  path  chosen 
in  the  2-plane  were  the  unit  circle  exactly. 

Of  course,  we  assume  that  this  image  is  actually  unknown  but  that  the  autocorrelation 
function  is  given  or  readily  calculated, 

4  2  13  6  3 

12  21  12  12  9 

r[m,  nj  =  7  19  49  19  7  (5.29) 

9  12  12  21  12 

3  6  13  2  4 

We  can  similarly  calculate  the  polynomial  associated  with  r[m.n],  pr(w.z).  The  paths  «»,■(<) 
in  the  w -plane  such  that  p,(wi(t).  z(t))  =  0  are  shown  in  Figure  5.13.  Note  that  the  paths 
shown  in  Figure  5.13  consist  of  those  corresponding  to  original  image,  drawn  as  a  solid  line, 
and  those  corresponding  to  the  image  rotated,  drawn  as  a  dashed  line.  The  present  algorithm 
does  not  calculate  all  these  paths  in  the  w-plane;  rather,  it  begins  with  a  one  zero  of  pf  (u>,  2(0)) 
and  generates  the  single  path  which  begins  at  this  zero.  Figure  5.14.  One  may  at  first  glance 
believe  that  there  might  be  some  ambiguity  at  the  location  in  the  tw-plane  where  two  paths 
cross.  However,  the  reader  should  note  that  an  ambiguity  will  only  occur  if  the  paths  cross  in 
the  ie-plane  at  the  same  value  of  z  and  w.  In  other  words,  one  may  visualize  these  paths  as 
trajectories  of  infinitesimal  particles  which  travel  over  the  tw-plane.  An  ambiguity  which  would 
preclude  the  possibility  of  isolating  a  single  path  would  only  occur  if  two  particles  collide,  i.e., 
are  on  the  same  position  at  the  same  time.  Using  samples  of  the  path  pair  (w(t).  z(t)).  the 
linear  equation  (5.22)  is  solved,  and  after  normalization,  results  in  the  image. 

1  0  3 

3  3  0  (5.30) 

1  2  4 

Comparing  Figures  5.11.  5.12  and  5.14.  we  see  that  we  followed  a  path  which  corresponds  to  the 
mirror  polynomial  of  pr(w,z ).  Thus,  the  resulting  reconstruction  is  rotated  by  180  degrees. 

The  next  section  is  concerned  with  some  issues  involved  in  implementing  the  algo¬ 
rithm  just  presented  for  images  much  larger  than  that  of  this  simple  example. 


5.3  Implementation  Issues 

The  theory  developed  in  the  previous  chapter  has  been  applied  to  write  a  computer 
program  to  extract  an  image  from  its  Fourier  transform  magnitude.  In  developing 
this  program  several  practical  considerations  need  to  be  addressed;  they  all  revolve,  as 
would  be  expected,  around  the  effects  of  finite  precision  representation  of  numbers  in 
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Figure  5.13:  Zero  paths  of  the  polynomial  associated  with  the  autocorrelation  function 
of  the  example.  Zeros  of  solution  polynomial  are  shown  as  solid,  those  corresponding 
to  the  mirror  polynomial  are  dashed. 
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a  computer.  In  the  section  below  we  discuss  these  considerations  and  how  they  have 
been  incorporated  into  the  final  algorithm. 

The  first  consideration  relates  to  the  proper  choice  for  a  path  z(t).  Recall  that 
in  the  description  of  the  algorithm,  z{t )  is  allowed  to  be  arbitrary,  as  long  as  it  is 
composed  strictly  of  ordinary  points  of  the  algebraic  function  defined  by  the  polynomial 
associated  with  the  autocorrelation  function,  pr(w,z).  Theoretically,  z(t)  could  be  a 
short  line  segment  in  the  z-piane.  In  this  case,  the  resulting  set  of  zero  samples  would 
be  representative  of  a  small  region  in  the  total  set  of  zeros  of  this  factor.  When  one  then 
tries  to  solve  the  set  of  linear  equations  given  by  (5.22)  to  extract  the  coefficients  of  this 
factor,  the  equation  coefficients  will  be  almost  the  same,  inviting  numerical  problems. 
The  objective  then  is  to  generate  a  path  z(t)  which  is  as  varied  as  possible  so  that  the 
resulting  samples  of  the  zeros  of  the  factor  polynomial  reflect  in  an  accurate  manner 
the  total  set  of  zeros  of  this  factor. 

A  second  consideration  concerns  the  evaluation  of  high  degree  polynomials.  In 
evaluating  pr[w,z),  z  has  to  be  raised  to  a  high  power,  namely  the  degree  of  pr{w,z) 
in  z.  For  an  N  by  N  image,  this  degree  is  2(iV  —  1).  If  the  magnitude  of  z  is  large 
or  small,  this  will  lead  to  severe  loss  in  numerical  precision.  This  consideration  then 
implies  that  the  magnitude  of  z(t)  for  all  t  should  be  in  the  proximity  of  one. 

Although  we  would  like  to  have  the  magnitude  of  z(t )  be  close  to  unity,  the  unit  circle 
itself  has  to  be  avoided.  Consider  a  signal  zjm,  nj  with  Fourier  transform  magnitude 
which  becomes  zero  at  (exp(ju),  exp{jv)).  Then  z  =  exp(ju)  and  w  =  exp{jv )  is  a 
zero  of  both,  p,( w,  z)  and  p,(tu,  z),  i.e.,  z  =  exp{ju)  is  a  singular  point.  Thus,  for  some 
images  there  are  singular  points  on  the  unit  circle  in  the  z-piane.  An  adequate  strategy 


is  to  pick  a  path  z(t)  where  the  magnitude  of  z(t)  is  some  number  ciose  to  but  not  equai 
to  one.  Once  z(t )  and  the  initial  value  w0  have  been  specified,  the  resulting  path  w(t) 
is  completely  determined.  For  this  reason,  it  cannot  be  assured  that  the  magnitude  of 
each  sample  w(t)  will  be  close  to  unity  even  if  the  magnitude  of  z(t)  is  close  to  unity. 
At  first  glance  this  seems  like  a  severe  problem,  and  in  fact  it  can  lead  to  numerical 
problems  in  some  cases.  However,  in  the  great  majority  of  cases  we  considered,  this 
problem  did  not  arise.  There  are  two  explanations  for  the  well-behavedness  of  u>(£); 
first,  w(t)  will  become  large  only  if  z(t)  passes  near  the  type  of  singular  point  where  the 
algebraic  function  t v{z)  becomes  infinite.  We  have  already  shown  that  there  are  few 
such  points.  Similarly,  one  can  show  that  w(t)  may  become  zero  only  if  z(t)  crosses  one 
of  a  finite  number  of  points  in  the  z-plane.  As  long  as  the  z(t)  path  does  not  come  close 
to  this  set  of  points,  we  can  expect  that  w(t)  will  be  well-behaved.  Second,  it  is  well 
known  that  the  roots  of  large  polynomials  with  arbitrary  coefficients  tend  to  cluster 
around  the  unit  circle  [41],  This  behavior  has  indeed  been  observed  in  the  numerical 
experiments. 

With  the  above  considerations  in  mind,  the  path  shown  in  Figure  5.15  was  chosen. 
The  path  z(t)  travels  counterclockwise  in  a  full  circle  of  radius  slightly  larger  than  1 
(ranging  from  1.1  to  1.2)  and  then  travel  clockwise  in  a  circle  of  radius  slightly  less 
than  1  (ranging  from  .8  to  .9).  At  this  point,  the  path  pair  is  “exchanged”.  Suppose 
that  at  the  end  of  the  depicted  path,  ine  zero  pair  is  (tu0,  -9).  At  this  point,  a  new 
path  is  begun  using  w  as  the  independent  variable.  This  path  begins  at  u>0,  goes  to 
w  =  1.1  and  then  travels  in  the  same  trajectory  as  z(t),  Figure  5.16.  As  long  as  the 
exchange  occurs  at  an  ordinary  point  of  z(w),  and  a  path  consisting  of  only  ordinary 
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points  is  used,  the  resulting  zeros  of  p,(ui,  z)  calculated  will  continue  to  correspond  to 
a  single  irreducible  factor  of  p,(tu,  z).  Intuitively,  this  “exchange”  between  w  and  z 
leads  to  an  equal  consideration  of  the  vertical  and  horizontal  orientations  of  the  image 
and  thus  presents  a  more  representative  sampling  of  the  zeros  of  one  factor  of  p*( w,  z) 
or  p,(w,z).  To  demonstrate  the  complexity  of  the  paths  involved,  in  Figures  5.17  and 
5.18  we  show  a  portion  of  a  typical  path  z(t)  and  the  resulting  path  t v(t),  for  a  20  by 
20  image. 

The  strategy  described  in  the  previous  paragraph  was  not  successful  in  ail  trials. 
We  have  applied  our  algorithm  to  about  40  images  of  sizes  20  by  20  and  larger,  and 
have  found  that  for  this  range  of  image  size  and  one  out  of  every  ten  images,  the  path 
picked  was  close  enough  to  a  singular  point  that  the  zeros  calculated  changed  from  one 
factor  to  another. 

One  way  to  avoid  this  is  to  calculate  where  all  the  singular  points  are  via  the 
discriminant  polynomial  discussed  in  Appendix  B,  and  then  generate  a  path  which  is 
far  from  any  singular  point.  An  alternate  strategy  is  to  keep  monitoring  the  partial 
derivative  of  pr(w,z)  with  respect  to  w.  If  this  derivative  becomes  zero,  then  the 
trajectory  may  be  approaching  a  singular  point.  Then,  the  location  of  the  singular 
point  can  be  located  and  circumvented.  We  have  not  implemented  either  of  these 
strategies  which  would  have  significantly  increased  the  computational  requirements  of 
the  program.  Since  this  situation  occurred  so  seldomly,  a  more  efficient  strategy  is  to 
restart  the  algorithm  using  a  slightly  different  path.  In  all  our  experiments  it  was  never 
found  necessary  to  use  more  than  two  trial  paths. 

The  next  consideration,  now  that  the  choice  of  the  path  z(t )  has  been  discussed,  is 


the  solution  of  the  differential  equation, 


dw{t)  _  pr,z(wt z)  dz(t) 
dt  pf'W(w,z)  dt 


(5.31) 


This  equation  is  converted  into  a  system  of  two  simultaneous  ordinary  differential 
equations  by  considering  the  real  and  imaginary  parts  separately.  Equation  (5.31)  then 
becomes  a  (nonlinear)  vector  initial  value  problem  with  a  governing  differential  equa¬ 
tion  of  first  order.  This  problem  has  been  studied  extensively  and  several  algorithms 
have  been  developed  which  use  variable  step  size  techniques  to  achieve  a  pre-specified 
accuracy  [42].  In  practice,  we  have  found  it  more  computationally  efficient  to  use  the 
differential  equation  solver  to  find  a  few  digits  of  accuracy  for  w  given  z  and  then  use 
the  following  version  of  Newton’s  method  to  calculate  the  zero  of  pr(w,z)  accurately. 
Suppose  that  an  estimate  of  a  zero  of  pt(w,z)  for  a  given  z  has  been  found,  w0.  This 


estimate  can  then  be  improved  via  the  iteration, 


Wi+l  =  Wi  - 


Pr{vi,z) 

Pr.A  Wi.*) 


(5.32) 


One  last  issue  to  be  resolved  is  the  proper  computation  of  one  solution  to  the 
homogeneous  set  of  equations  (5.22).  Recall  that  the  associated  matrix  ,  say  M,  has 
rank  deficiency  of  1,  i.e.,  it  has  a  unique  (to  a  constant)  non-trivial  null  vector.  The 
approach  we  have  taken  is  to  use  the  QR  algorithm  with  column  pivoting  [43].  This 
algorithm  calculates  the  following  factorization  of  an  m  by  n  <  m  matrix  M  of  rank 


n  —  1: 


M  =  QRU 


(5.33) 


.  v.  •- 


where  Q  is  an  m  by  m  orthogonal  matrix,  R  is  of  the  form 


Ri,i  t  n  -  1 

0  0  1 


(5.34) 


n  -  1  1 


where  RlA  is  non-singular  upper  triangular,  r  is  a  column  vector,  0  is  a  row  vector 
consisting  of  all  zeros,  and  IT  is  a  permutation  matrix.  Once  we  have  the  factorization 
above,  we  can  find  a  solution  to  Mi  =  0  as 


y  =  -Rulr 


(5.35) 


i  =  n 


(5.36) 


Chapter  6 

Numerical  Results 


In  this  chapter  we  illustrate  the  performance  of  the  new  phase  retrieval  algorithm 
by  applying  it  to  several  families  of  images.  We  have  applied  this  algorithm  to  about 
40  different  realistic  images  of  size  20  by  20  and  above  but  we  will  just  be  showing 
some  representative  examples. 

We  follow  these  examples  by  a  study  of  the  effect  of  noisy  observations  on  the 
reconstruction  algorithm  and  a  comparison  of  the  computational  requirements  of  the 
new  algorithm  with  the  Deighton  algorithm. 

6.1  Examples  of  Application  of  Factorization  Algo¬ 
rithm  to  Phase  Retrieval 

In  this  section  we  present  several  examples  of  the  application  of  the  phase  re¬ 
trieval  algorithm  to  several  different  types  of  images.  The  experimental  procedure  is  as 
follows.  The  autocorrelation  function  of  the  original  image  is  first  calculated  and  the 
polynomial  associated  with  this  autocorrelation  function  formed.  The  phase  retrieval 
algorithm  is  then  used  to  factor  this  polynomial.  The  coefficients  of  the  irreducible  fac¬ 
tor  extracted  then  correspond  to  the  pixel  values  of  the  reconstructed  image.  In  every 


case  it  is  assumed  that  the  support  of  the  image,  or  rather,  the  size  of  the  smallest  rect¬ 
angle  surrounding  the  support  of  the  original  image,  is  known.  In  order  to  emphasize 
the  fact  that  we  are  actually  reconstructing  an  image  from  its  autocorrelation  function, 
we  display  the  original  image,  its  autocorrelation  function,  and  the  reconstruction.  In 
all  examples,  it  is  assumed  that  the  unknown  image  corresponds  to  a  set  of  intensi¬ 
ties  which  must  be  non-negative.  This  requirement  is  certainly  not  necessary  for  our 
algorithm  to  be  applicable;  it  only  removes  the  sign  ambiguity. 

In  Figure  6.1  we  applied  the  phase  retrieval  algorithm  to  reconstruct  a  25  by  25  pixel 
image.  This  size  image  is  the  largest  that  could  be  accommodated  in  our  Vax/1 1-750 
computer  system  storage  capability  without  excessive  memory  paging.  In  this  case, 
the  reconstruction  has  the  same  orientation  as  the  original.  In  the  second  example, 
another  25  by  25  pixel  image,  Figure  6.2  was  reconstructed.  In  this  example,  the 
resulting  reconstruction  is  rotated  by  180  degrees.  As  mentioned  earlier  since  for  the 
case  of  rectangular  support,  an  image  and  its  180  degree  rotated  version  have  the  same 
autocorrelation  function,  it  is  impossible  to  restore  the  absolute  orientation  of  the 
image.  Figures  6.3  and  6.4  show  the  result  of  applying  the  reconstruction  algorithm 
to  different  families  of  square  images.  Figure  6.3  is  a  high  contrast  image  of  size  20 
by  20  pixels  which  may  be  representative  of  an  application  in  astronomy.  Note  that  in 
this  case,  the  support  of  the  image  is  not  assumed  known,  since  it  is  not  rectangular; 
however,  it  is  assumed  that  the  smallest  rectangle  surrounding  the  support  is  known. 
The  image  in  Figure  6.4  consists  of  a  20  by  20  image  composed  of  independent  random 
noise  uniformly  distributed  between  0  and  1.  This  example  is  especially  interesting 
since  the  autocorrelation  function  of  images  with  random  coefficients  are  very  similar. 


Figure  6.1:  Girl  image,  (a)  Original  Image,  (b)  Autocorrelation  function,  (c)  Recon¬ 


struction. 
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This  is  due  to  the  fact  that  most  of  the  autocorrelation  function  content  is  concentrated 
at  r[0,  Oj.  This  is  also  the  family  of  objects  for  which  the  iterative  algorithms  performed 
very  poorly. 

If  the  support  of  the  image  is  non-symmetric  so  that  the  image  rotated  by  180 
degrees  would  not  “fit”  in  the  known  support  constraint,  then  the  rotation  ambiguity  is 
removed.  However,  the  algorithm  described  in  this  thesis  does  not  use  this  information. 
Implicit  in  our  development  is  that  extracting  either  p2[xv,  2)  or  p,(u/,  z)  is  adequate. 
In  the  situation  of  a  triangular  support,  p,(u;,  z)  is  associated  with  a  signal  which  does 
not  satisfy  the  support  constraint.  However,  our  algorithm  is  just  as  likely  to  return 
such  an  invalid  signal.  For  example,  in  Figure  6.5  we  show  the  result  of  applying  our 
algorithm  to  a  20  by  20  pixel  image  with  triangular  support.  The  reconstruction  (by 
chance)  has  the  wrong  support,  since  it  is  rotated  by  180  degrees.  This  effect  however, 
is  of  minor  concern  since  it  is  a  trivial  task  to  extract  the  correct  reconstruction  from 
such  a  rotated  version.  This  is  a  second  example  where  only  a  bound  on  the  support 
is  assumed  known,  rather  than  the  exact  bound. 

6.2  Effect  of  Noise  on  Reconstruction 

In  an  effort  to  characterize  the  susceptibility  of  the  algorithm  to  errors  in  the  known 
data,  we  studied  the  effect  of  adding  independent  gaussian  noise  to  the  known  autocor¬ 
relation  function  of  a  real  image,  and  then  applying  the  phase  retrieval  algorithm  to  this 
noisy  autocorrelation  function.  Adding  noise  in  this  domain  is  equivalent  in  some  sense 
to  adding  noise  to  the  Fourier  transform  magnitude  squared.  Since  most  applications 
of  phase  retrieval  measure  either  the  autocorrelation  function  or  the  Fourier  transform 
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magnitude  squared,  it  was  felt  that  this  form  of  degradation  is  more  pertinent  than 
adding  noise  to  the  Fourier  transform  magnitude  directly. 

Consider  the  autocorrelation  function  r[m,  n]  of  an  unknown  image  x [m,  n].  A  noisy 
version  of  r[m,n]  was  calculated,  f[m,  n]  via, 

f[m,n]  =  r[m,n]  +  tu[m,n)  (6.1) 


The  signal  i v{m,  nj  consisted  of  simulated  gaussian  random  noise  of  variance  <x?  with 
the  added  condition  that  w[-m,  -nj  =  tn[m,  nj.  Thus,  it  is  assumed  that  only  non- 
redundant  measurements  of  the  autocorrelation  function  were  made. 

The  resulting  autocorrelation  function  f[m,  nj  is  thus  symmetric  about  the  origin. 
Thus,  this  degradation  assumes  that  the  signal  is  known  to  be  real  and  furthermore, 
that  the  support  of  the  image,  i.e.,  the  support  of  the  true  autocorrelation  function, 
is  known  exactly.  Given  a  correlation  signal  to  noise  value,  SNR< ,  the  variance  of  the 
added  noise  was  calculated  as 


2  _  Em.n  r[m,n]2 
‘  SNR 


(6.2) 


The  phase  retrieval  algorithm  was  then  applied  using  this  noisy  autocorrelation  func¬ 
tion.  The  reconstruction  signal  to  noise  ratio,  SNR0,  was  calculated  as, 

SNR0  =  ~  *lm’nD2 

£m.nzjm,nj2 

where  x\m,n)  is  the  original  image  and  z[m,  nj  is  the  reconstructed  image.  The  SNRn 
was  maximized  over  a  rotation  of  the  reconstruction  by  180  degrees. 


The  experiments  were  performed  for  two  image  sizes,  20  by  20,  and  16  by  16.  For 
each  image  size,  two  images  were  used,  denoted  here  as  “1”  and  “2”,  and  two  sets  of 


pseudo-random  noise  fields,  “A”  and  “BB.  The  magnitude  of  the  noise  was  changed 
until  SNR0  defined  in  (6.3)  was  in  the  range  of  10.  The  results  of  these  experiments 


are  displayed  in  Figures  6.6  and  6.7. 


Reconstruction  SNR 


Data  SNR 

Figure  6.6:  Effect  of  noise  on  reconstruction  for  16  by  16  pixel  images.  Graphs  shown 
correspond  tc  the  application  of  noise  signals  “A”  and  “B”  to  images  “1"  and “2”. 


We  note  that  the  algorithm  is  very  sensitive  to  noise.  It  is  difficult  to  specify 
the  sensitivity  of  the  algorithm  since  it  varies  dramatically  between  images  and  noise 
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signals,  (see  for  example  1-A  and  1-B  in  Figure  6.6).  However,  it  seems  that  for  16  by 
16  images,  a  S N R<  of  better  than  101#  is  required  for  faithful  reconstruction.  For  20 
by  20  images,  a  5JVJ2,  of  1019  seems  to  be  required  for  most  images. 

The  sensitivity  of  the  algorithm  to  small  amounts  of  noise  can  be  explained  with 
the  help  of  the  discussion  in  Section  3.2.1  regarding  reducible  and  irreducible  polyno¬ 
mials.  A  small  perturbation  to  the  true  associated  polynomial  pr(w,z)  will  result  in 
an  irreducible  polynomial.  As  the  perturbation  increases,  pr(w,z)  will  no  longer  be 
even  “approximately”  reducible.  This  will  result  in  difficulty  in  finding  a  good  approx¬ 
imation  to  the  factors  of  the  original  factorable  polynomial  pr(w,z).  A  more  formal 
argument  supporting  this  observation  is  presented  in  [44]. 

6.3  Computational  Requirements 

In  order  to  compare  the  computational  requirements  of  the  new  algorithm  with 
Deighton’s  algorithm,  we  measured  the  number  of  CPU  seconds  which  are  required 
to  reconstruct  a  typical  image  using  the  new  algorithm.  The  solid  plot  in  Figure 
6.8  is  a  display  of  the  number  of  CPU  seconds  required  per  reconstruction  plotted 
against  image  size.  We  see  that  for  small  images  the  Deighton  algorithm  is  much 
more  efficient.  However,  for  image  sizes  larger  than  about  11,  the  Deighton  algorithm 
becomes  computationally  more  expensive. 

We  have  found  experimentally  that  as  the  image  size  increases,  the  QR  decomposi¬ 
tion  becomes  a  dominant  part  of  the  computation  in  the  new  algorithm.  This  leads  us 
to  hypothesize  that  the  algorithm  asymptotically  requires  computations  in  the  order  of 
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jV6  for  reconstructing  an  N  by  N  image.  In  contrast,  the  Deighton  algorithm  requires 
computations  which  grow  exponentially  with  N.  A  more  efficient  implementation  of 
either  the  Deighton  or  our  algorithm  would  certainly  affect  the  point  at  which  our 
algorithm  performs  better  than  the  Deighton  algorithm.  However,  the  important  thing 
to  note  is  the  rate  of  growth  of  each  algorithm  rather  than  the  actual  computer  usage. 

We  can  also  compare  the  storage  requirements  of  the  Deighton  algorithm  against 
the  present  algorithm.  For  the  Deighton  algorithm,  the  lists  which  need  to  be  compared 
have  sizes  which  depend  exponentially  on  image  size.  In  the  case  of  the  new  algorithm, 
the  storage  requirements  are  set  predominantly  by  the  set  of  linear  equations  to  be 
solved.  The  storage  requirements  for  this  matrix  are  on  the  order  of  2 N*  for  an  N 
by  N  image.  As  the  image  size  increases,  the  storage  requirements  for  the  Deighton 
algorithm  increase  much  faster  than  for  the  new  algorithm. 


CPU  seconds 

200G0  - 


o  new  algorithm 


15000 


o  deighton  algorithm 

10000 


o 


5000  • 


Image  size 


Figure  6.8:  Comparison  of  CPU  usage  for  Deighton  vs.  new  algorithm. 


Chapter  7 


Other  Applications 


There  are  other  applications  which  may  benefit  from  the  ideas  developed  in  this  the¬ 
sis.  We  discuss  here  two  such  situations.  The  first  involves  the  general  problem  of  fac¬ 
torization  of  polynomials  in  two  variables,  while  the  second  addresses  two-dimensional 
recursive  filter  stability  testing. 


7.1  Application  to  General  Bivariate  Polynomial 
Factorization 

As  an  example  of  where  such  an  algorithm  may  be  applicable,  consider  a  two- 
dimensional  linear  shift-invariant  system  which  is  specified  by  the  difference  equation 

[451, 

H  H  buy'™  -  k,  n  -  1}  =  H  dkix[m  -  k,  n  -  1}  (7.1) 


where  y[m,  n]  is  the  system  response  to  an  excitation  x[m,nj.  The  z-transform  of  the 
response  of  this  system  to  a  unit  impulse  can  be  expressed  as, 


„(  v  _  LkEt  aktw-kz-‘  _  A(w,z) 
{  ’  J  E* Ei bkiw-kz~l  B(w,z) 


(7.2) 


We  see  that  A{w,  z )  and  B(w,  z )  can  be  expressed  as  two  bivariate  polynomials  in  the 
variables  w~l,  z~l.  In  characterizing  H{ w,  z)  one  has  to  answer  the  question  of  whether 
the  representation  given  in  (7.2)  is  minimal,  i.e.,  whether  there  are  some  polynomials 
A'{w,z )  and  B'(w,z)  of  total  degree  less  than  A(w,z)  and  B(w,  z)  respectively  which 


satisfy, 


H(w,z) 


A'{w,z) 
B'(w,  z) 


This  is  equivalent  to  asking  whether  A(w,  z )  and  B{ w,  z)  are  such  that 


(7.3) 


A(w ,  z)  =  P(w,  z)A'(w,  z )  (7.4) 

B(w,  z)  =  P[w,  z)B'(w,  z)  (7.5) 


where  P(w,z)  is  a  non-constant  polynomial.  In  other  words,  one  may  ask  whether 
A(w,z)  and  B(w,z)  are  relatively  prime  (coprime).  Several  algorithms  have  been 
developed  to  test  whether  two  polynomials  are  relatively  prime  [46,47],  However,  these 
tests  are  very  complicated  for  large  polynomials,  since  effectively,  they  generate  the 
resultant  R.\d{z )  described  in  Appendix  B. 

A  necessary  condition  for  the  required  P{w,z)  to  exist  is  that  A(w,z)  and  B{w,z ) 
be  factorable.  Thus,  by  using  a  polynomial  bivariate  factorization  algorithm,  we  can 
check  whether  coprimeness  is  satisfied.  If  either  A(tn,  z)  or  B(w,z)  is  irreducible,  then 
there  is  no  need  to  use  the  more  complicated  coprimeness  algorithms. 
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There  are  some  modifications  required  from  the  aigorithm  presented  in  Chapter  5. 
These  modifications  are  due  to  the  fact  that,  unlike  the  phase  retrieval  problem,  one 
does  not  know  a  priori  the  degree  of  the  factor  isolated  by  the  zero-tracking  procedure. 
Thus,  although  one  can  set  a  bound  on  the  maximum  number  of  homogeneous  equations 
which  need  to  be  solved,  one  does  not  know  a  priori  the  number  of  unknowns,  i.e.,  the 
number  of  coefficients  in  the  polynomial  to  be  isolated.  In  Appendix  C,  we  discuss  two 
methods  for  estimating  the  total  degree  of  the  isolated  factor.  Here  we  demonstrate 
the  use  of  the  second  algorithm  discussed  in  Appendix  C. 

To  illustrate  the  performance  of  our  factorization  algorithm,  we  have  generated  a 
test  polynomial  p(w,z )  given  by, 


z  i  2  8  13  23  30  12 

4  13  18  38  48  17 

3  12  23  48  55  16 

3  14  25  44  40  10 

2  8  17  25  19  4 

1  5  7  10  6  1 


which  can  be  expressed  as  the  product  of  the  following  two  polynomials, 


z  i  2  2  I 

6  5  4 

3  2  1 


In  this  example,  the  total  degree  of  p(w,  z)  is  10.  From  Bezout’s  theorem,  we 
know  that  any  other  polynomial  of  total  degree  less  than  or  equal  to  10  which  is  prime 
relative  to  p(w,  z )  can  have  at  most  100  zeros  in  common  with  p(w,  z ).  Thus  we  need  to 
calculate  at  most  101  zeros  along  a  path  pair  ( w(t),z(t ))  consisting  of  ordinary  points. 
This  way  we  are  assured  that  these  zeros  correspond  to  a  single  irreducible  factor  of 
p(w,z).  The  following  matrix  is  then  generated, 


1 

zi 

■  •  z10 
Z1 

ZiUl! 

■  ZyW* 

•  tu}° 

M  = 

1 

*2 

Z2 

••  z4° 

* 

Z2U>2 

z2wl 

(101  by  66)  (7.9) 

1 

z101 

zfoi  ' 

■  •  zin 
Z101 

ZlOl^lOl 

■■  *101^01  • 

•  • 

wi0l 

Note  that  although  the  degree  of  p(w,  z )  is  known,  we  need  to  assume  a  triangular 
support  for  the  coefficients  of  p(w,  z).  This  is  because  we  can  only  bound  the  total 
degree  of  the  factor  of  p(tu,  z),  not  necessarily  its  degree  in  each  variable.  Thus  although 
p(u!,z)  has  36  non-zero  coefficients,  M  has  66  columns.  According  to  Appendix  C,  M 
will  have  a  rank  deficiency,  i.e.,  66  -  r(Af),  equal  to 

(10  -  totdeg(d)  +  1)(  10  -  totdeg(d)  -I-  2)/2  (7.10) 


where  d(w,z)  is  the  irreducible  factor  of  p(w,z)  isolated  by  our  procedure.  In  order 
to  calculate  the  rank  of  M,  we  performed  a  QR  decomposition  of  M;  Table  7.1  shows 
some  of  the  diagonal  elements  of  the  triangular  matrix  R. 


diagonal  element 

value 

48 

0.00534025 

49 

0.00507684 

50 

-0.00145097 

51 

-0.00059735 

52 

4.45087e-15 

53 

-1.2814e-14 

54 

8.82197e-15 

55 

5.67471e-15 

.  _  _  i 

Table  7.1:  Diagonal  elements  of  triangular  matrix  after  QR  decomposition  for  calcu¬ 
lating  total  degree  of  factor. 

From  this  table  we  see  that  the  rank  of  M  is  51.  Therefore,  the  total  degree  of 
the  factor  d(w,z)  isolated  from  p{w,z )  is  6.  Once  the  total  degree  of  d(w,z)  has  been 
determined,  we  can  regenerate  M  with  fewer  columns, 


1 

*1 

o 

■  z 6 
*1 

ZitVi 

•  •  ZyXlll 

xu\ 

M  = 

1 

Z2 

z\ 

•  Z2 

Z2W2 

■  Z2ui\ 

■  •  w\ 

(101  by  15) 

1 

Zioi 

zioi 

■■  zm 

*101«>101 

••  2j01^?01  • 

••  <01 

(7.11) 

A  new  QR  decomposition  of  this  matrix  yields  a  single  non-trivial  null  vector  (to  a 
constant)  which  can  then  be  associated  with  the  coefficients  of  d(w,z),  given  below, 


w 

l  -.25 

-.25 

2.6e  - 

16 

-.25  — 2.3e  —  17  -8.3e  -  17 

-.25 

2.1e  -  16 

-.25 

-.25  5.9e  —  16  1.3e  -  16 

-.5 

-.5 

-.5 

-.5  -9.7c -17 

-1 

- 

—  .  i  0 

-  5 

-.25 

-2.1e  -  16 

5.7e  -  16 

4.0e  - 

16 

—  2.8e  -  16 

2.9e  -  16 

— 2.1e  -  16 

(7.12) 


We  see  that  to  a  scale  constant,  the  retrieved  factor  is  equal  to  one  of  the  generating 


polynomials. 


The  second  example  below  shows  how  the  algorithm  can  be  used  to  test  irreducibil- 
ity. 

For  this  example  we  consider  the  polynomial, 

w 

z  i  1  0  10  1 

01010  (7.13) 

10  10  1 
0  10  10 
10  10  1 

In  this  case  the  total  degree  of  this  polynomial  is  8,  so  a  65  by  45  matrix  M  of  the 
form  given  by  (7.9)  is  generated.  Table  7.2  shows  the  last  5  diagonal  elements  of  the 
triangular  matrix  calculated  via  a  QR  decomposition  of  M  for  this  case. 


diagonal  element 

value 

41 

-0.119653 

42 

-0.11235 

43 

0.110724 

44 

0.0845307 

45 

4.17397e-13 

Table  7.2:  Diagonal  elements  of  triangular  matrix  after  QR  decomposition  for  irre- 
ducibility  test. 


Since  the  rank  deficiency  of  M  is  seen  to  be  1  in  this  case,  we  conclude  that  the 
polynomial  shown  in  (7.13)  is  irreducible. 


7.2  Application  of  Root-tracking  to  Two-Dimensional 
Recursive  Filter  Stability  Testing 


A  further  application  of  the  ideas  presented  in  this  thesis  lies  in  the  field  of 


u 


stability  theory  of  two-dimensional  recursive  filters.  It  is  well  known  [45]  that  the 


question  of  whether  a  filter  with  rational  z-transform  transfer  function, 


H(tv,z)  = 


A(w,z) 

B(xvtz) 


(7.14) 


is  stable  or  not  depends  on  the  location  of  the  zeros  of  the  denominator  polynomial 
B(w,z)  l.  Specifically,  it  has  been  shown  that  if  all  branches  of  tw(z)  for  z  =  exp(jv), 
called  root  maps  in  the  literature,  are  of  magnitude  less  than  one,  then  the  filter 
described  in  (7.14)  is  bounded-input,  bounded-output  stable.  The  most  difficult  part 
of  this  test  is  showing  that  the  branches  of  w(z)  do  not  cross  the  unit  circle  in  the 
iu-plane,  i.e.,  that  A(exp(ju),  exp(jv))  is  never  zero. 


A  computationally  attractive  algorithm  for  testing  for  filter  stability  relies  on  the 
fact  that  the  two-dimensional  unwrapped  phase  must  be  continuous  and  periodic  if 
all  root  maps  are  always  less  than  one  in  magnitude  [48].  This  procedure  is  fast  and 
efficient,  as  long  as  all  the  root  maps  are  far  from  the  unit  circle.  However,  as  the  root 
maps  get  close  to  the  unit  circle,  the  phase  becomes  almost  discontinuous,  rendering  a 
phase  unwrapping  algorithm  helpless. 

At  this  point  it  becomes  more  convenient  to  look  at  the  root  maps  themselves, 
rather  than  consider  their  effect  on  the  phase.  Shaw  [48]  proposes  detecting  when  the 
unwrapped  phase  is  almost  discontinuous  and  then  doing  a  constrained  search  for  a 
minimum  of  the  magnitude  of  A{exp(ju),  exp(jv)).  However,  for  complicated  filters, 
the  magnitude  function  is  likely  to  have  local  minima  leading  to  an  erroneous  answer  or 


'Strictly  speaking,  it  is  necessary  to  consider  nonessential  singularities  of  the  second  kind  where 
;)  and  D(w.  ;)  are  zero  simultaneously  on  the  unit  bicircle.  However,  this  is  extremely  rare  and  is 
usually  iguored  in  stability  testing  procedures. 
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almost  flat  regions,  which  results  in  extremely  slow  convergence  of  the  search  algorithm. 
In  fact,  in  Shaw’s  example,  the  search  for  a  local  minimum  took  14  times  longer  than 
the  phase  unwrapping  itself. 

As  an  alternative,  we  propose  that  the  root  maps  be  tracked  directly  for  those 
frequency  regions  where  the  phase  is  almost  discontinuous.  In  other  words,  calculating 
accurately  the  branch  of  w(z)  for  z  =  exp( v0)  to  2  =  exp(vi)  in  the  range  of  v,  v0  <  v  < 
Vi  for  which  the  phase  is  almost  discontinuous.  Such  a  procedure  would  provide  exact 
knowledge  of  whether  the  filter  is  stable,  and  if  it  is  stable,  how  close  are  the  zeros  to 
the  unit  circle,  thus  providing  a  margin  of  stability. 

The  actual  calculation  of  the  root  maps  for  all  values  of  2  =  exp(j  v)  is  also  useful 
as  an  aid  for  testing  stability  by  itself.  It  provides  a  way  of  graphically  portraying  the 
margin  of  stability  of  a  filter,  especially  for  small  filters,  such  as  8  by  8. 

As  an  illustration  of  the  possible  application  of  this  algorithm  to  filter  stability 
testing,  we  used  the  differential  equation  method  of  Chapter  5  to  calculate  all  the  root 
maps  of  the  following  filter, 


w 

z  l  1.000000  -1.11870  .355  (7.15) 

-1.000030  1.55704  -.550496 

.400020  -.659803  .222204 

The  polynomial  corresponding  to  A(tu,  1)  was  calculated  and  its  roots  determined. 
Each  of  these  roots  provided  an  initial  condition  for  the  differential  equation  solver 
Samples  of  the  root  map  were  thus  calculated  and  are  displayed  in  Figure  7.1. 

From  the  samples  of  the  root  maps  calculated  it  was  inferred  that  a  root  map  may 
cross  the  unit  circle.  The  root  map  was  recalculated  in  a  smaller  interval  centered 
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around  the  region  where  the  root  map  may  cross  the  unit  circle.  The  resulting  root 
map  is  shown  in  Figure  7.2.  Since  this  precise  root  map  does  not  cross  the  unit  circle, 

§{tt/} 
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0.8- 

0.7- 

0.6  — 

0.5  — 

0.5  0.6  0.7  0.8  0.9  • 

u;-plane 

Figure  7.2:  Single  root  map  of  two-dimensional  recursive  filter  over  small  region. 

we  can  be  confident  that  the  filter  is  stable. 


Chapter  8 


Summary  and  Suggestions  for 
Future  Research 

The  object  of  this  thesis  has  been  to  develop  a  better  algorithm  for  reconstruction 
of  two-dimensional  discrete  signals  from  the  Fourier  transform  magnitude.  To  this  effect 
we  have  developed  a  new  closed-form  algorithm  for  reconstruction  which  is  guaranteed 
to  yield  a  correct  solution  given  accurate  data  and  which  is  computationally  much  more 
efficient  than  previous  closed-form  algorithms,  both  in  terms  of  time  of  execution  and 
memory  requirements.  The  algorithm  has  been  applied  to  a  large  number  of  images 
successfully. 

We  have  noted  however,  that  memory  requirements  preclude  its  use  for  images  larger 
than  25  by  25  in  the  computer  facilities  available  to  us.  Furthermore,  the  algorithm  is 
extremely  sensitive  to  noise  in  the  given  Fourier  transform  magnitude  information. 

Many  algorithms  have  been  proposed  in  the  mathematical  literature  for  solving 
large  sets  of  linear  equations  in  a  memory  efficient  manner;  for  example,  the  conjugate 
gradients  method.  Such  an  approach  may  be  fruitful  in  reducing  the  memory  require¬ 
ments  of  our  phase  retrieval  procedure.  This  would  allow  the  reconstruction  of  much 


larger  images. 


The  problem  of  noise  sensitivity  is  an  important  one.  It  is  not  clear  at  this  point 
whether  the  sensitivity  problem  is  particular  to  our  algorithm  only  or  is  shared  by 
other  closed-form  algorithms.  In  a  sense,  this  property  is  characteristic  of  the  phase 
retrieval  problem  since  once  noise  is  added  to  the  Fourier  transform  magnitude  or 
the  autocorrelation  signal,  the  reconstruction  problem  is  replaced  by  an  approximation 
problem,  where  no  solution  will  satisfy  the  data  exactly.  In  this  case  we  need  to  decide 
on  a  criterion  of  goodness  of  match  between  a  candidate  solution  and  the  measured 
data.  Such  a  criterion  should  be  amenable  to  analysis  and  the  development  of  an 
appropriate  closed-form  algorithm.  The  GSF  algorithm  of  Fienup  does  reduce  an  error 
norm;  however,  it  may  not  drive  the  norm  to  its  minimum  value.  In  our  opinion,  solving 
the  phase  retrieval  approximation  problem  should  be  a  major  focus  of  future  research. 

In  developing  our  phase  retrieval  algorithm  we  have  considered  it  as  strictly  a  poly¬ 
nomial  factorization  problem.  We  made  little  use  of  the  fact  that  the  factors  of  pr(ti/,  z), 
Pj(w,z )  and  z),  are  not  independent  but  are  intricately  related.  The  use  of  this 

symmetry  may  prove  useful  in  advancing  the  ideas  presented  here. 

Finally,  we  have  presented  a  few  possible  applications  of  the  new  factorization  al¬ 
gorithm  and  root-tracking  procedure  to  the  area  of  two-dimensional  filter  design.  It 
would  be  very  interesting  to  develop  other  applications  of  these  ideas. 


Appendix  A 

Convergence  of  Iterative 
Algorithms  for  Phase  Retrieval 


In  this  appendix  we  conduct  a  study  of  the  convergence  properties  of  the  Fienup 
algorithms  for  reconstruction  from  Fourier  transform  magnitude.  To  achieve  this  aim 
we  have  conducted  a  Monte  Carlo  study,  i.e.,  the  application  of  each  algorithm  to  a 
large  number  of  randomly  generated  images.  The  desire  for  a  large  sample  size  dictated 
that  a  small  image  size  be  used.  However,  such  a  study  does  provide  a  qualitative  idea 
of  the  behavior  of  the  algorithms  in  general. 

Each  trial  consisted  of  the  generation  of  two  random  images;  the  first  image  was 
used  as  the  target  image  to  be  reconstructed,  while  the  second  was  used  as  an  initial 
estimate.  There  were  four  sets  of  experiments,  studying  the  effect  of  a  symmetric  vs. 
non-symmetric  support  constraint  and  image  size.  The  parameters  of  each  of  the  four 
experiments  is  summarized  in  Table  A.l. 


In  the  first  two  experiments,  “A”  and  “B”.  the  images  consisted  of  a  4  by  4  field 
of  independently  generated  random  numbers  uniformly  distributed  between  0  and  1. 
It  has  been  noted  that  the  iterative  algorithms  perform  better  when  a  nonsymmetric 


experiment 

image  size 

support 

max.  iterations 

times  repeated 

4x4 

rect. 

600 

100 

4x4 

triang. 

600 

100 

C 

8x8 

rect. 

400 

50 

D 

8x8 

triang. 

400 

50 

Table  A.l:  Summary  of  parameters  used  in  convergence  study  of  iterative  phase  re¬ 
trieval  algorithms. 


region  of  support  is  known,  for  example,  if  the  support  is  triangular.  To  study  the 
effect  of  a  nonsymmetric  support,  experiment  “A”  used  a  rectangular  support,  while 
experiment  “B”  used  a  triangular  support  by  setting  the  the  lower  triangular  half  of  the 
image  to  zero.  The  target  image  and  the  initial  estimate  were  generated  independently 
and  no  effort  was  made  to  generate  a  good  initial  estimate,  other  than  satisfying  the 
known  spatial  domain  constraints.  Thus,  we  are  only  interested  in  the  case  where 
the  only  a  priori  information  is  the  Fourier  transform  magnitude,  the  signal  support, 
and  possibly,  that  the  signal  is  non-negative.  From  the  target  image,  zero  padded  to 
an  8  by  8  image,  the  discrete  Fourier  transform  was  calculated  and  the  magnitude 
used  as  information  to  the  reconstruction  algorithms.  Each  algorithm  was  run  with 
the  same  target  and  initial  estimate  images.  For  each  algorithm,  an  error  measure 
was  calculated  between  the  estimated  image  at  the  midpoint  iteration  and  the  target 
image.  The  algorithm  was  then  continued  and  the  error  measure  was  recalculated  at 
the  final  iteration.  The  trials  were  repeated,  with  different  target  and  initial  estimates, 
100  times.  The  error  measure  used  was  the  normalized  mean  squared  error  (NMSE) 


between  the  estimate  z,[m,  nj  and  the  target  image  z[m,  nj. 


NMbh  r.v-i  r.V-u.im 


(A.l) 


The  NMSE  was  minimized  over  all  trivial  associates  of  z,[m,  n],  i.e.,  the  estimate  was 


successively  multiplied  by  -1,  and/or  rotated  by  180  degrees. 

In  order  to  study  the  effect  of  image  size  on  the  quality  of  reconstruction,  the 
experiments  described  above  were  repeated  for  images  consisting  of  a  rectangular  or 
triangular  field  of  size  8  by  8,  resulting  in  experiments  “C”  and  “D”.  In  this  case,  a 
total  of  400  iterations  were  used  and  the  experiment  was  repeated  50  times.  Again  the 
NMSE  was  recorded  both  at  the  midpoint  and  final  iteration. 

After  conducting  these  experiments  we  were  able  to  discern  certain  properties  of  the 
iterative  algorithms  studied.  Calculating  the  NMSE  at  the  middle  iteration  as  well  as 
the  final  iteration  allows  us  to  discern  if  the  algorithm  has  progressed  substantially  in 
the  intervening  iterations.  A  clear  way  to  compare  the  performance  of  these  algorithms 
at  the  midpoint  and  final  iteration  is  to  draw  a  scatter  plot  where  the  abscissa  of  each 
point  is  the  NMSE  value  at  the  midpoint  iteration  while  the  ordinate  is  the  NMSE 
value  at  the  final  iteration.  Points  which  lie  close  to  the  main  diagonal  indicate  that 
no  change  has  occurred  in  the  intervening  iterations  while  points  substantially  below 
the  main  diagonal  indicate  substantial  improvement  between  the  midpoint  and  final 
iterations. 

Figures  A.l  to  A. 3  show  scatter  plots  of  the  NMSE  at  the  midpoint  and  finai 
iteration  for  each  algorithm  during  experiment  “A”.  A  similar  set  of  scatter  plots  are 
displayed  in  Figures  A. 4  to  A. 6  for  experiment  “B”. 

In  the  above  figures,  we  see  that  for  NMSE  greater  than  10— 4 ,  there  is  no  change 
between  the  midpoint  and  final  iteration  in  almost  all  trials.  Therefore,  we  can  assume 
that  for  large  NMSE  (  greater  than  10-4)  the  algorithms  have  reached  a  fixed  point  of 
the  iterative  procedure  by  the  final  iteration,  and  thus  cannot  be  expected  to  make  any 


600th  iter 


Figure  A.l:  Error  at  600th  vs.  300th  iteration  for  GSF  -  rectangular  4  by  4.  Number 
in  box  denotes  number  of  trials  which  fell  in  that  square. 
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Figure  A. 3:  Error  at  600th  vs.  300tb  iteration  for  HIO  -  rectangular  4  by  4 
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Figure  A. 4:  Error  at  600th  vs.  300th  iteration  for  GSF  -  triangular  4  by  4 


133 


V, 

\ *. 


600tn  iter 


D.  i 


O.Oi 


0-001 


0.000' 


i*-05 


i®-06 


ift-07 


'e-00 


’e*09 


V  w  V  "®  w— ■  -5— w  'W  —  T,  -» 


B* 


3 

-4 

>3 

1 

>1 

1 

J 

i 

j 

'  J 

1 


►  -  .- 


[3I 


i  51 


300t*  Hr 


>-'3  -e-’?  i®-m  ^e-’Q  ie-09  ’c-08  ie-97  ie-06  ie-05  0.0001  0.00i  0.01  0. 

e*-ro'-  a*  SCO*.*'  ’te-  vs.  330t^  iter  for  GSf-p  -  triangular  4  t>y  4 


h\*’ 

Ky 


Figure  A. 5:  Error  at  600th  vs.  300th  iteration  for  GSF-P  -  triangular  4  by  4 
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Figure  A. 6:  Error  at  600th  vs.  300th  iteration  for  HIO  -  triangular  4  by  4 


more  progress  in  subsequent  iterations.  The  above  observations  are  consistent  with  the 
results  of  previous  authors  ( 19,2,2 lj  who  have  noted  the  slow  progress  of  the  GSF  and 
GSF-P  algorithms. 

Confidence  that  a  fixed  point  of  the  iterative  procedure  has  been  reached  allows 
us  to  evaluate  the  success  rate  of  each  algorithm.  For  this  purpose,  we  need  to  use  a 
criterion  for  success.  Although  an  appropriate  cutoff  value  for  an  acceptable  NMSE 
is  difficult  to  assess,  especially  for  small  images  as  those  considered  in  this  paper,  we 
have  used  a  cutoff  of  10'2  and  10"4,  which  correspond  to  approximately  one  and  two 
significant  figures  for  each  pixel  value,  respectively.  For  this  cutoff,  the  success  rate  of 
each  algorithm  for  4  by  4  images  is  noted  in  Table  A.2. 


[  Success  rate  of  convergence  (percent)  | 

support  type 

rectangular 

triangular  | 

21  9 

59 

53 

22  12 

60 

55 

HIO 

22  12 

60 

52 

Table  A. 2:  Summary  of  algorithm  convergence  as  a  function  of  support  and  criterion 
for  4  by  4  image. 

The  experiments  described  above  were  repeated  for  8  by  8  images,  experiments 
“C”  and  “D”.  Scatter  plots  of  the  results  are  shown  in  Figures  A. 7  to  A.  12.  These 
plots  allow  us  to  again  conclude  that  for  the  vast  majority  of  trials,  a  fixed  point  of 
the  iteration  has  been  reached  by  the  final  iteration.  The  success  rate  of  the  iterative 
algorithms  for  8  by  8  images  is  summarized  in  Table  A.3. 

Conclusions  derived  from  Tables  A. 2  and  A.3  are  discussed  in  the  main  text. 
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Figure  A.9:  Error  at  400tb  vs.  200th  iteration  for  HIO  -  rectangular  8  by  8 
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Figure  A.  10:  Error  at  400th  vs.  200th  iteration  for  GSF  -  triangular  8  by  8 
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Figure  A.  12:  Error  at  400th  vs.  200th  iteration  for  HIO  -  triangular  8  by  8 
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Appendix  B 

A  Bound  on  the  Number  of  Finite 
Common  Zeros  Based  on 
Polynomial  Degree  1 


Recall  that  Bezout’s  theorem  is  concerned  with  determining  the  number  of  com¬ 
mon  zeros  of  two  bivariate  polynomials.  It  is  restated  here  for  convenience, 

Theorem  B.l  If  p(w,z)  and  q(w,z)  are  bivariate  polynomials  of  total  degree  r  and 
s  respectively  with  no  common  factors,  then  there  are  at  most  rs  distinct  pairs  (w,z) 
where 

p(w,z)  =  q{w,z)  =  0  (B.l) 

Since  Bezout’s  theorem  is  concerned  with  total  degree  as  opposed  to  degree  in  each 
variable,  it  pertains  most  generally  to  polynomials  whose  coefficients  have  triangular 
support  as  shown  in  Figure  B.l. 

In  our  case  of  reconstruction  from  Fourier  transform  magnitude,  this  corresponds 
to  an  image  which  has  a  triangular  support.  On  the  other  hand,  one  is  many  times 
interested  in  images  with  square  or  rectangular  support  as  shown  in  Figure  B.2. 

For  the  case  when  the  polynomials  under  consideration  have  rectangular  support, 

'The  derivation  presented  here  is  due  to  a  collaboration  between  A.  Zakhor  and  the  author  and  is 
also  described  in  [49!. 
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for  a  polynomial  of  degree  (2,2). 


we  are  able  to  lower  the  bound  on  the  number  of  common  finite  zeros  from  the  bound 
set  by  Bezout’s  theorem.  Specifically  if  the  two  relatively  prime  polynomials  p(w,z) 
and  q(w,  z)  are  given  by 

A/p  .Vp 

PK*)  =  E  E  Pm,numzn  (B.2) 

m~0  n=0 

and 

A/,  -V, 

=  E  E  (B-3) 

m=0  n=0 

the  upper  bound  on  the  number  of  common  finite  zeros  set  by  Bezout’s  theorem  is 
(A',  +  Mp){Nq  +  Mq).  Our  objective  is  to  establish  a  tighter  upper  bound  on  the 
number  of  common  finite  zeros  of  p{w,z)  and  q(w,z). 

Before  proceeding,  we  need  to  review  several  results  concerning  the  resultant  of  poly¬ 
nomials  in  one  or  two  variables.  The  resultant  Rpq  of  two  one-dimensional  polynomials 
p(u;)  and  q(w) 

Mp 

(B.4) 


A/p 

pM  =  E  P"U;" 

n=0 

A/, 

=  E  9"W" 

n=0 


(B.5) 

(B.6) 


is  defined  [50]  as  the  determinant  of  the  [Mp  +  M,)  by  ( Mp  +  Mq)  matrix 
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A  basic  property  of  resultants  is  stated  in  the  following  theorem  [50] , 


Theorem  B.2  When  the  polynomials  p(tu)  and  q(w)  have  numerical  coefficients,  a 
necessary  and  sufficient  condition  that  they  shall  have  a  finite  or  infinite  common  root 
is  that  =  0. 


Consider  now  the  two  relatively  prime  bivariate  polynomials  p{w,z)  and  q{w,z) 
defined  in  (B.2)  and  (B.3)  expressed  as  polynomials  in  u>  with  coefficients  which  are 
polynomials  in  z , 


p{w,z)  =  po(z)-t- Pi (*)«>  +  •• 

•  +  Pm,{z)vF’ 

(B.8) 

<7(10.  z)  =  q0[z)  +  <?i(z)tu  + 

+  9m,  (*)«'*'* 

(B.9) 

We  can  define  the  resultant  of  p(tu,  z )  and  q(w,  z)  with  respect  to  in  as  the  determinant 
of  the  (Mp  +  Mq)  by  (Mp  4-  Mq)  matrix,  M(z),  with  polynomial  entries: 
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(B.  10) 

This  resultant  is  a  function  of  the  remaining  variable  z  and  is  denoted  by  Rpq{z)- 


Expanding  the  determinant  of  the  above  matrix,  and  taking  into  account  that  each 
p,(z)  and  q,(z )  is  of  degree  at  most  Np  and  Nq  respectively,  we  can  conclude  that 
Rpq{z)  is  a  polynomial  of  degree  NpMq  4-  NqMP  or  less  It  can  be  shown  ibj,  that  if 
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p[w,  z )  and  q(w,  z )  are  relatively  prime  then  Rpq{z)  is  not  identically  zero.  Moreover, 
if  (wo,  Zo)  is  a  common  zero  of  p(w,  z)  and  q(w,  z)  then  Rpq(zo)  =  0.  Thus  the  zero  sets 
of  p(vj,  z)  and  q(w,  z)  have  at  most  NpMq  +  NqMp  values  of  z  in  common. 

As  our  argument  stands,  we  have  not  yet  placed  any  tight  limit  on  the  number  of 
intersections  of  p(w,  z )  and  q(w,  z )  since  for  each  z,  which  is  a  root  of  Rpq{z )  there  could 
be  a  large  number  of  wt,  such  that  for  each  j, 

p{w})z,)  =  q(w},z,)  =  0  (B.  11) 

In  order  to  specify  the  number  of  Wj  for  each  z^,  we  need  to  study  the  behavior  of 
iZ,l(J(z)  in  the  vicinity  of  each  z,. 

Theorem  B.3  If  at  each  z0  there  are  k  values  of  tv,  w},  such  that 

P{v>j,z o)  =  o)  =  0  for  j  =  1,  •  •  • ,  k  (B.12) 

then  Rf,q{z)  has  a  zero  of  multiplicity  k  at  z0. 

The  above  theorem  implies  that  p( w,  z )  and  z)  as  defined  in  equation  (B.2) 
and  (B.3)  have  at  most  NpMq  +  NqMp  zeros  in  common. 

In  order  to  prove  Theorem  B.3  we  need  to  review  some  results  on  matrices  with 
polynomial  entries,  relating  to  the  Smith  Normal  Form  [51], 

Theorem  B.4  Let  A{x)  be  an  n  by  n  polynomial  matrix  of  normal  rank  r.  We  can 
find  untmodular  matrices  {P[x),Q{x)) ,  such  that 

B(z)  =  P(x)A{x)Q{z)  (B.13) 

and 

1.  B(x)  is  a  diagonal  polynomial  matrix  called  the  Smith  Normal  Form  of  A(x). 
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2.  The  first  r  diagonal  elements  of  B(x)  are  monte  polynomials  Pi(x),  p2(x),  ■, 

Pr(z)- 

S.  The  remaining  diagonal  elements,  if  any,  are  zero. 

4.  pi( x)  divides  p<+ i(x)  for  i  =  1,  •  •  • ,  r  -  1. 

Normal  rank  in  the  above  theorem  is  defined  in  [51].  The  unimodular  polynomial 
matrices  of  the  above  theorem  are  defined  to  have  nonzero  constant  determinant  inde¬ 
pendent  of  x.  Therefore,  if  r  =  n,  i.e.,  if  A{ x)  has  full  normal  rank  then, 

det(B(x))  =  JJ  Pi(x)  (B.  14) 

1=1 

Also,  from  part  (4)  of  Theorem  B.4  we  can  conclude  that  if  p,(z)  =  0  then  p*(x)  =  0 
for  k  >  i.  From  the  above  theorem,  we  can  derive  the  following, 

Theorem  B.5  Let  A(x)  be  a  polynomial  matrix  of  full  normal  rank.  If  A(z0)  has  rank 
deficiency  of  k  then  det(A(x)),  has  a  zero  of  multiplicity  k  at  x0. 

Proof:  Using  Theorem  B.4  we  can  find  B(x),  the  Smith  normal  form  of  A(x).  Since 
P{x)  and  Q(x)  are  unimodular,  B(x)  is  of  full  normal  rank.  Furthermore,  the  ranks  of 
B(x)  and  A{x)  at  each  value  of  x,  including  x0,  are  equal.  Therefore  B(x0)  has  rank 
deficiency  of  k.  This  means  that  p„(z 0)  =  Pn-i(io)  =  •  ■  •  =  pn-t+i(*o)  =  0.  Therefore, 
from  (B.14)  det(B(z))  has  a  zero  of  multiplicity  k  at  x0.  Since  the  determinant  of  A{x) 
is  within  a  constant  factor  of  that  of  J9(x),  A{x)  also  has  a  zero  of  multiplicity  k  at  z0. 
C 

Using  Theorem  B.5,  we  can  return  to  Theorem  B.3, 

Proof:  ( Theorem  B  S)  Suppose  that  for  z0  there  are  k  common  finite  zeros  to,,  j  = 
1  ...  A;  between  p(uj,z)  and  q(w,  z).  Then  the  matrix  M(x0)  defined  by  (B.  10)  must 


for  ;  =  1,  • ,  A:.  Thus  M(z0)  has  rank  deficiency  of  k.  Furthermore,  since  p(w,  z )  and 

q(w,  z)  have  no  common  factors,  Rpq(z)  is  not  identically  zero,  so  M(z)  has  full  normal 
rank.  From  Theorem  B.5  then,  Rpq(z)  has  a  zero  of  multiplicity  k  at  z0.  □ 

From  Theorem  B.3,  and  the  fact  that  Rpq(z)  is  a  polynomial  of  degree  NpMq+NqMp, 
we  get  immediately,  the  main  result  of  this  appendix, 


Theorem  B.6  Let  p(u>,z)  and  q(w,z )  be  two  polynomials  of  degree  ( MP,NP )  and 
( Mq,Nq )  with  no  common  factors,  then  there  are  at  most  NpMq  +  NqMp  pairs  of  finite 
complex  numbers  ( w ,  z )  such  that 

p(w,z)  =  q(w,z)  =0  (B-16) 
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Appendix  C 

Calculating  the  Degree  of  an 
Irreducible  Factor 

In  this  appendix  we  want  to  consider  the  situation  where  the  degree  of  an  irre¬ 
ducible  factor  of  the  polynomial  p(w,  z )  is  not  known.  From  the  path-tracking  algorithm 
developed  in  Chapter  5,  we  have  generated  a  large  number  of  zeros  of  the  polynomial 
p(w,z).  We  denote  these  zeros  by  (u>*,Zi)  for  k  =  1,  -,K.  We  will  see  later  that 

for  this  case  K  should  be  picked  to  be  the  maximum  allowable  by  Bezout’s  theorem, 
totdeg{p)2  +  1.  From  our  earlier  arguments,  these  zeros  all  correspond  to  a  single  ir¬ 
reducible  factor  of  p(w,z),  which  we  will  denote  by  d(w,  z).  It  is  this  factor  which  we 
seek  to  isolate.  Although  the  total  degree  of  p(w,  z)  is  known,  we  only  know  that  the 
total  degree  of  d(w,  z)  can  be  at  most  totdeg(p).  In  this  appendix  we  describe  two 
mechanisms  for  calculating  totdeg(d). 

The  first  mechanism  is  the  simplest  conceptually.  Assume  that  totdeg(d)  =  1.  Then 
for  all  zeros  ( ),  there  must  be  constants  doo,dOi,dl0  which  are  not  all  zero  such 
that 

doo  +  doiU>k  +  dl0Zif  =  0  (C.l) 
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for  all  k.  Equivalently,  we  need  that  the  matrix 

1  W\  Z\ 

1  XV2  z2 

•  •  • 

•  •  • 

1  WK  ZK 

have  a  non-trivial  nullspace,  i.e.,  a  non-zero  vector  i  such  that  Afii  =  0.  Thus,  the 

first  step  of  this  procedure  is  to  calculate  the  rank  of  Mu  which  we  will  symbolize  by 

r(M,).  If  it  is  less  than  3,  then  totdeg(d)  =  1.  On  the  other  hand,  if  has  full  rank, 

then  we  generate  M2  defined  by, 

- 

1  U)!  xv\  Zi  tWjZ!  z\ 

1  w2  u>l  z2  w2z2  z\ 

M2=  (C.3) 

•  •  • 

•  •  • 

•  •  • 

1  wK  w2K  Zk  WkZk  z2k 

Again,  the  rank  of  M2  is  calculated.  If  it  is  less  than  6,  then  totdeg{d)  ~  2,  otherwise, 
we  calculate  M3.  Thus  the  algorithm  proceeds  by  successively  calculating  the  rank  of 
A/„  and  checking  if  the  rank  of  M„  is  less  than  (n  -I-  l)(n  +  2)/2.  The  smallest  n  for 
which  this  condition  is  true  is  the  total  degree  of  d(w,z). 

For  large  polynomials,  the  procedure  described  above  becomes  very  inefficient  in 
some  cases.  For  example,  if  the  polynomial  p(w,  z)  is  itself  irreducible,  totdcg(p)  differ¬ 
ent  matrices  and  rank  determinations  would  be  required. 

It  turns  out  that  by  examining  the  rank  of  only  the  last  matrix  M„  for  n  =  totdeg{p) 
we  can  calculate  the  degree  of  d{w,  z).  For  simplicity,  we  will  call  this  Mn  simply  M  for 
the  rest  of  the  discussion.  Since  by  assumption  d[wi,,Zk)  -  0  for  all  k,  the  coefficients 


of  d(w,  z),  with  appropriately  padded  zero  coefficients,  form  a  null  vector  of  M.  In  fact, 
the  coefficients  of  any  polynomial  of  total  degree  totdeg(p)  which  has  d(ty,  2)  as  a  factor 
form  a  null  vector  of  M.  Conversely,  with  any  null  vector  of  A/,  we  can  associate  a 
“null”  polynomial  with  d{ w,  2)  as  a  factor  by  virtue  of  Bezout’s  theorem.  The  problem 
now  is  to  find  the  dimension  of  this  nullspace.  Consider  the  set  of  polynomials, 

d[w,  2),  zd[w,  2),  wd{w,  2),  •  •  ■ ,  wmznd(w,  z)  (C.4) 

for  all  m,  n  such  that  the  total  degree  of  wmznd(w,  z)  is  less  than  or  equal  to  totdeg(p). 
Then  all  the  null  polynomials  can  be  expressed  as  a  linear  combination  of  these  polyno¬ 
mials.  In  fact  the  “coordinates”  of  a  null  polynomial  q(w,  2 )  =  d[w,  z)a(tv,  2)  are  exactly 
the  coefficients  of  a(w,  2).  Thus  the  set  of  polynomials  in  (C.4)  span  the  nullspace  of  M. 
We  have  not  shown  however,  that  these  polynomials  are  linearly  independent.  Suppose 
that  they  form  a  linearly  dependent  set,  then  we  can  generate  a  set  of  coefficients  c.,y 
not  all  zero  such  that 

Y  c,  ,xu'zid{w,  2)  =  0  (C.5) 

«■/ 

In  other  words,  there  exists  a  polynomial  c(w,z)  such  that  c{w,z)d{w,z)  =  0  for  all 
10,2.  However,  since  d(w,  2)  is  not  identically  zero,  c(xv,  z)  must  be  zero.  Thus  the 
polynomials  in  (C.4)  form  a  linearly  independent  set  which  spans  the  nullspace  of  M. 
From  this  argument  we  see  that  the  dimension  of  the  nullspace  is  simply  the  number 
of  polynomials  in  the  set  described  in  (C.4),  which  in  turn  is  simply  the  maximum 
number  of  coefficients  in  the  polynomial  c(u/,z)  in  (C.5).  The  total  degree  of  c(u;,z)  is 
totdeg(p)  -  totdeg(d).  A  polynomial  of  total  degree  n  has  (n  4-  l)(n  -I-  2)/2  coefficients; 


hence,  the  dimension  of  the  nullspace  of  M  is 


( totdeg(p )  —  totdeg(d)  4-  1  )(totdeg(p)  -  totdeg(d)  +  2)/2  (C.6) 


Note  that  not  all  possible  values  for  the  rank  of  M  are  possible,  since  some  values 
for  the  rank  of  M  lead  to  non-integer  values  of  totdeg(d).  The  procedure  to  calculate 
totdeg(d)  can  be  summarized  as  follows: 

1.  Generate  M  consisting  of  ( totdeg(p )  4-  1  )(totdeg(p)  4-  2)/2  columns  and  K  = 
totdeg(p)2  4-  1  rows. 


2.  Calculate  the  rank  of  M,  r(M). 

3.  The  value  of  totdeg(d)  is  given  by 

,  .  .  2 totdeg(p)  +3  4-  yj{2 totdeg(p)  4-  3)2  -  4r(M) 

totdeg{d)  - - - - - -  (C.7) 

M 


Bibliography 


[lj  Yu.  M.  Bruck  and  L.  G.  Sodin.  On  the  ambiguity  of  the  image  reconstruction 
problem.  Optics  Communications ,  30:304-308,  1979. 

[2]  M.  H.  Hayes.  The  reconstruction  of  a  multidimensional  sequence  from  the  phase 
or  magnitude  of  the  Fourier  transform.  IEEE  Trans.  Acoustics,  Speech,  and  Signal 
Processing,  ASSP-30:140-154,  1982. 

[3]  J.  R.  Fienup.  Space  object  imaging  through  the  turbulent  atmosphere.  Optical 
Engineering,  18:529-534,  1979. 

[4]  A.  Labeyrie.  Attainment  of  diffraction  limited  resolution  in  large  telescopes  by 
Fourier  speckle  patterns  in  star  images.  Astron.  and  Astrophys.,  6:85-87,  1970. 

[5]  R.  W.  Gerchberg  and  W.  O.  Saxton.  Phase  determination  from  image  and  diffrac¬ 
tion  plane  pictures  in  the  electron  microscope.  Optik,  34:275-284,  1971. 

[6]  G.  N.  Ramachandran  and  R.  Srinivasan.  Fourier  Methods  in  Crystallography. 
Wiley-Interscience,  New  York,  1970. 

[7]  A.  V.  Oppenheim  and  R.  W.  Schafer.  Digital  Signal  Processing.  Prentice  Hall, 
Englewood  Cliffs,  NJ,  1975. 

[8j  A.  Mostowski  and  M.  Stark.  Introduction  to  Higher  Algebra.  Pergamon  Press, 
New  York,  1964. 

[9]  M.  H.  Hayes.  Signal  Reconstruction  from  Phase  or  Magnitude.  PhD  thesis,  Mas¬ 
sachusetts  Institute  of  Technology,  1981. 

[lOj  G.  A.  Bliss.  Algebraic  Functions.  American  Mathematical  Society,  New  York, 
1933. 

[llj  E.  J.  Akutowicz.  On  the  determination  of  the  phase  of  a  Fourier  integral.  Trans. 
American  Mathematical  Society,  83:179-182,  1956. 

[12]  A.  Walther.  The  question  of  phase  retrieval  in  optics.  Optica  Acta,  10:41-49,  1963. 

[13]  E.  M.  Hofstetter.  Construction  of  time-limited  functions  with  specified  autocorre¬ 
lation  functions.  IEEE  Trans.  Information  Theory,  IT-8:1 19-126,  1964. 


[14]  M.  H.  Hayes  and  J.  H.  McClellan.  Reducible  polynomials  in  more  than  one  vari¬ 
able.  Proceedings  of  the  IEEE,  70:197-198,  1982. 

[15]  J.  L.  C.  Sanz,  T.  S.  Huang,  and  F.  Cukierman.  Stability  of  unique  Fourier  trans¬ 
form  phase  reconstruction.  Journal  of  the  Optical  Society  of  America,  73:1442- 
1445,  1983. 

[16]  A.  M.  J.  Huiser  and  P.  van  Toorn.  Ambiguity  of  the  phase-reconstruction  problem. 
Journal  of  the  Optical  Society  of  America,  46:407-420,  1976. 

[17]  J.  L.  C.  Sanz  and  T.  S.  Huang.  Unique  reconstruction  of  a  band-limited  multidi¬ 
mensional  signal  from  its  phase  or  magnitude.  Journal  of  the  Optical  Society  of 
America,  73:1446-1450,  1983. 

[18]  J.  R.  Fienup.  Reconstruction  of  an  object  from  the  modulus  of  its  Fourier  trans¬ 
form.  Optics  Letters,  3:27-29,  1978. 

[19]  J.  R.  Fienup.  Phase  retrieval  algorithms:  A  comparison.  Applied  Optics,  21:2758- 
2769,  1982. 

[20]  R.  W.  Gerchberg  and  W.  O.  Saxton.  A  practical  algorithm  for  the  determination 
of  phase  from  image  and  diffraction  plane  pictures.  Optik,  35:237-246,  1972. 

[21]  J.  L.  C.  Sanz,  T.  S.  Huang,  and  T-.  F.  Wu.  A  note  on  iterative  Fourier  transform 
phase  reconstruction  from  magnitude.  IEEE  Trans.  Acoustics,  Speech,  and  Signal 
Processing,  ASSP-32:1251-1254,  1984. 

[22]  A.  Levi  and  H.  Stark.  Image  restoration  by  the  method  of  generalized  projections 
with  application  to  restoration  from  magnitude.  Journal  of  the  Optical  Society  of 
America  -  A,  1:932-943,  1984. 

[23]  M.  C.  Won,  D.  Mnyama,  and  R.  H.  T.  Bates.  Improving  initial  phase  estimates 
for  phase  retrieval  algorithms.  Optica  Acta,  32:377-396,  1985. 

[24]  J.  R.  Fienup.  Comments  on  ‘The  reconstruction  of  a  multidimensional  sequence 
from  the  phase  or  magnitude  of  the  Fourier  transform’.  IEEE  Trans.  Acoustics, 
Speech,  and  Signal  Processing,  ASSP-31:1256-1262,  1983. 

[25]  R.  H.  T.  Bates.  Fourier  phase  problems  are  uniquely  solvable  in  more  than  one 
dimension.  I.  Underlying  theory.  Optik,  61:247-262,  1982. 

[26]  K.  L.  Garden  and  R.  H.  T.  Bates.  Fourier  phase  problems  are  uniquely  solvable  in 
more  than  one  dimension.  II.  One-dimensional  considerations.  Optik,  62:131-142, 
1982. 

[27]  W.  R.  Fright  and  R.  H.  T.  Bates.  Fourier  phase  problems  are  uniquely  solvable  in 
more  than  one  dimension.  III.  Computational  examples  for  two  dimensions.  Optik, 
62:219-230.  1982. 


[28]  R.  H.  T.  Bates  and  W.  R.  Fright.  Composite  two-dimensionai  phase-restoration 
procedure.  Journal  of  the  Optical  Society  of  America,  73:358-365,  1983. 

[29]  N.  Canterakis.  Magnitude-only  reconstruction  of  two-dimensional  sequences  with 
finite  regions  of  support.  IEEE  Trans.  Acoustics,  Speech,  and  Signal  Processing, 
ASSP-31:1256-1262,  1983. 

[30]  H.  V.  Deighton,  M.  S.  Scivier,  and  M.  A.  Fiddy.  Solution  of  the  two-dimensional 
phase- retrieval  problem.  Optics  Letters,  10:250-251,  1985. 

[31]  H.  M.  Berenyi,  H.  V.  Deighton,  and  M.  A.  Fiddy.  The  use  of  bivariate  polynomial 
factorization  algorithms  in  two-dimensional  phase  problems.  Optica  Acta,  32:689- 
701,  1985. 

[32]  M.  Marden.  Geometry  of  Polynomials.  American  Mathematical  Society,  Provi¬ 
dence,  Rhode  Island,  1966. 

[33]  A.  I.  Markushevich.  Theory  of  Functions  of  a  Complex  Variable.  Volume  II, 
Prentice-Hall,  Englewood  Cliffs,  NJ,  1965. 

[34]  E.  Kaltofen.  Factorization  of  polynomials.  Computing  Supplement ,  4:95-113, 
1982. 

[35]  E.  Kaltofen.  A  polynomial-time  reduction  from  bivariate  to  univariate  integral 
polynomial  factorization.  In  Proc.  2Srd  Symposium  on  Foundations  of  Computer 
Science,  IEEE,  pages  57-64,  1982. 

[36]  J.  Lipson.  Newton’s  method:  a  great  algebraic  algorithm.  Proc  1976  Symposium 
on  Symbolic  and  Algebraic  Computation,  260-270,  1982. 

[37]  D.  Y.  Y.  Yun.  Hensel  meets  Newton  -  algebraic  constructions  in  an  analytic 
setting.  In  J.  F.  Traub,  editor,  Algebraic  Computational  Complexity,  Academic 
Press,  New  York,  1976. 

[38]  P.  J.  Napier  and  R.  H.  T.  Bates.  Inferring  phase  information  from  modulus  in¬ 
formation  in  two-dimensional  aperture  synthesis.  Astron.  Astrophys.  Supplement, 
15:427-430,  1974. 

[39]  S.  R.  Curtis,  J.  S.  Lim,  and  A.  V.  Oppenheim.  Signal  reconstruction  from  one  bit 
of  Fourier  transform  phase.  IEEE  Trans.  Acoustics,  Speech,  and  Signal  Processing, 
ASSP-33:643-657,  1985. 

[40]  S.  R.  Curtis.  Reconstruction  of  Multidimensional  Signals  from  Zero  Crossings. 
PhD  thesis,  Massachusetts  Institute  of  Technology,  1985. 

[41]  K.  Steiglitz  and  B.  Dickinson.  Phase  unwrapping  by  factorization.  IEEE  Trans. 
Acoustics,  Speech,  and  Signal  Processing,  ASSP-30:984-991,  1982. 

42]  G.  Dahlquist  and  A.  Bjorck.  Numerical  Methods.  Prentice-Hall,  Englewood  Cliffs, 
NJ,  1974. 


[43]  G.  H.  Golub  and  C.  F.  Van  Loan.  Matrix  Computations.  Johns  Hopkius  University 
Press,  Baltimore,  1983. 

[44]  J.  L.  C.  Sanz  and  T.  S.  Huang.  Polynomial  system  of  equations  and  its  application 
to  the  study  of  the  effect  of  noise  on  multidimensional  Fourier  transform  phase 
retrieval  from  magnitude.  IEEE  Trans.  Acoustics,  Speech,  and  Signal  Processing, 
ASSP-33:997-1004,  1985. 

[45]  D.  E.  Dudgeon  and  R.  M.  Mersereau.  Multidimensional  Digital  Signal  Processing. 
Prentice-Hall,  Englewood  Cliffs,  NJ,  1984. 

[46]  N.  K.  Bose.  Applied  Multidimensional  System  Theory.  Van  Nostrand  Reinhold, 
New  York,  1982. 

[47]  E.  Emre  and  O.  Huseyi  j.  Relative  primeness  of  multivariate  polynomials.  IEEE 
Trans.  Circuits  and  Systems ,  22:56-57,  1975. 

[48]  G.  A.  Shaw.  Design,  stability  and  performance  of  two-dimensional  recursive  digital 
filters.  PhD  thesis,  Georgia  Institute  of  Technology,  1979. 

[49]  A.  Zakhor  and  D.  Izraelevitz.  A  note  on  the  sampling  of  zero-crossings  of  two- 
dimensional  signals.  Accepted  for  publication,  March  1986. 

[50]  R.  J.  Walker.  Algebraic  Curves.  Dover  Publications,  New  York,  1962. 

[51]  T.  Kailath.  Linear  Systems.  Prentice-Hall,  Englewood  Cliffs,  NJ,  1980. 


mm  i  '  k  '  J  ^  *  ■  ■  ►  '  J  J  JA  j*  ja  J*  /  _  m  .  »  m 


DISTRIBUTION  LIST 


DODAAD 


Director  HX1241 

Defense  Advanced  Research  Project  Agency 

1400  Wilson  Boulevard 

Arlington,  Virginia  22209 

Attn:  Program  Management 


Head  Mathematical  Sciences  Division  N00014 

Office  of  Naval  Research 
800  North  Quincy  Street 
Arlington,  Virginia  22217 


Administrative  Contracting  Officer  N66017 

E19-628 

Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts  02139 


Director  N0017  3 

Naval  Research  Laboratory 
Attn:  Code  2627 
Washington,  D.  C.  20375 


Defense  Technical  Information  Center  S47  031 

Bldg  5,  Cameron  Station 
Alexandria,  Virginia  22314 


Dr.  Judith  Daly 
DARPA  /  TTO 
1400  Wilson  Boulevard 
Arlington,  Virginia  22209 


Code 


(1) 


(1) 


(1) 


(6) 


(12) 


(1) 


